home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / DCLAP 6d / dclap6d / SeqPups / appsrc / fastdnaml.src / fastDNAml.c < prev    next >
Text File  |  1996-07-05  |  147KB  |  5,203 lines
  1. #define programName        "fastDNAml"
  2. #define programVersion     "1.1.1a"
  3. #define programVersionInt   10101
  4. #define programDate        "November 2, 1994"
  5. #define programDateInt      19941102
  6.  
  7.  
  8.  
  9. /*  fastDNAml
  10.  *
  11.  *  Olsen, G. J., Matsuda, H., Hagstrom, R., and Overbeek, R.  1994.
  12.  *  fastDNAml:  A tool for construction of phylogenetic trees of DNA
  13.  *  sequences using maximum likelihood.  Comput. Appl. Biosci. 10: 41-48.
  14.  *
  15.  *
  16.  *  Based in (large) part on the program dnaml by Joseph Felsenstein.
  17.  *
  18.  *  Felsenstein, J.  1981.  Evolutionary trees from DNA sequences:
  19.  *  A maximum likelihood approach.  J. Mol. Evol. 17: 368-376.
  20.  *
  21.  *
  22.  *  Copyright notice from dnaml:
  23.  *
  24.  *  version 3.3. (c) Copyright 1986, 1990 by the University of Washington
  25.  *  and Joseph Felsenstein.  Written by Joseph Felsenstein.  Permission is
  26.  *  granted to copy and use this program provided no fee is charged for it
  27.  *  and provided that this copyright notice is not removed.
  28.  */
  29.  
  30. /*  Conversion to C and changes in sequential code by Gary Olsen, 1991-1994
  31.  *
  32.  *  p4 version by Hideo Matsuda and Ross Overbeek, 1991-1993
  33.  */
  34.  
  35. /*
  36.  *  1.0    March 14, 1992
  37.  *         Initial "release" version
  38.  *
  39.  *  1.0.1  March 18, 1992
  40.  *         Add ntaxa to tree comments
  41.  *         Set minimum branch length on reading tree
  42.  *         Add blanks around operators in treeString (for prolog parsing)
  43.  *         Add program version to treeString comments
  44.  *
  45.  *  1.0.2  April 6, 1992
  46.  *         Improved option line diagnostics
  47.  *         Improved auxiliary line diagnostics
  48.  *         Removed some trailing blanks from output
  49.  *
  50.  *  1.0.3  April 6, 1992
  51.  *         Checkpoint trees that do not need any optimization
  52.  *         Print restart tree likelihood before optimizing
  53.  *         Fix treefile option so that it really toggles
  54.  *
  55.  *  1.0.4  July 13, 1992
  56.  *         Add support for tree fact (instead of true Newick tree) in
  57.  *            processTreeCom, treeReadLen, str_processTreeCom and
  58.  *            str_treeReadLen
  59.  *         Use bit operations in randum
  60.  *         Correct error in bootstrap mask used with weighting mask
  61.  *
  62.  *  1.0.5  August 22, 1992
  63.  *         Fix reading of underscore as first nonblank character in name
  64.  *         Add strchr and strstr functions to source code
  65.  *         Add output treefile name to message "Tree also written ..."
  66.  *
  67.  *  1.0.6  November 20, 1992
  68.  *         Change (! nsites) test in setupTopol to (nsites == 0) for MIPS R4000
  69.  *         Add vectorizing compiler directives for CRAY
  70.  *         Include updates and corrections to parallel code from H. Matsuda
  71.  *
  72.  *  1.0.7  March 25, 1993
  73.  *         Remove translation of underlines in taxon names
  74.  *
  75.  *  1.0.8  April 30, 1993
  76.  *         Remove version number from fastDNAml.h file name
  77.  *
  78.  *  1.0.9  August 12, 1993
  79.  *         Version was never released.
  80.  *         Redefine treefile formats and default:
  81.  *             0  None
  82.  *             1  Newick
  83.  *             2  Prolog
  84.  *             3  PHYLIP (Default)
  85.  *         Remove quote marks and comment from PHYLIP treefile format.
  86.  *
  87.  *  1.1.0  September 3-5, 1993
  88.  *         Arrays of size maxpatterns moved from stack to heap (mallocs) in
  89.  *            evaluateslope, makenewz, and cmpBestTrees.
  90.  *         Correct [maxsites] to [maxpatterns] in temporary array definitions
  91.  *            in Vectorize code of newview and evaluate.  (These should also
  92.  *            get converted to use malloc() at some point.)
  93.  *         Change randum to use 12 bit segments, not 6.  Change its seed
  94.  *            parameter to long *.
  95.  *         Remove the code that took the absolute value of random seeds.
  96.  *         Correct integer divide artifact in setting default transition/
  97.  *            transversion parameter values.
  98.  *         When transition/transversion ratio is "reset", change to small
  99.  *            value, not the program default.
  100.  *         Report the "reset" transition/transversion ratio in the output.
  101.  *         Move global data into analdef, rawDNA, and crunchedDNA structures.
  102.  *         Change names of routines white and digit to whitechar and digitchar.
  103.  *         Convert y[] to yType, which is normally char, but is int if the
  104.  *            Vectorize flag is set.
  105.  *         Split option line reading out of getoptions routine.
  106.  *
  107.  *  1.1.1  September 30, 1994
  108.  *         Incorporate changes made in 1.0.A (Feb. 11, 1994):
  109.  *            Remove processing of quotation marks within comments.
  110.  *            Break label finding into copy to string and find tip.
  111.  *            Generalize tree reading to read trees when names are and are not
  112.  *               already known.
  113.  *            Remove absolute value from randum seed reading.
  114.  *         Include integer version number and program date.
  115.  *         Remove maxsite, maxpatterns and maxsp limitations.
  116.  *         Incorporate code for retaining multiple trees.
  117.  *         Activate code for Hasegawa & Kishino test of tree differences.
  118.  *         Make quick add the default, with Q turning it off.
  119.  *         Make emperical frequencies the option with F turning it off.
  120.  *         Allow a residue frequency option line anywhere in the options.
  121.  *         Increase string length passed to treeString (should be length
  122.  *               checked, but ...)
  123.  *         Introduce (Sept.30) and fix (Oct. 26) bug in bootstrap sampling.
  124.  *         Fix error when user frequencies are last line and start with F.
  125.  */
  126.  
  127. #ifdef Master
  128. #  undef  Master
  129. #  define Master     1
  130. #  define Slave      0
  131. #  define Sequential 0
  132. #else
  133. #  ifdef Slave
  134. #    undef Slave
  135. #    define Master     0
  136. #    define Slave      1
  137. #    define Sequential 0
  138. #  else
  139. #    ifdef Sequential
  140. #      undef Sequential
  141. #    endif
  142. #    define Master     0
  143. #    define Slave      0
  144. #    define Sequential 1
  145. #  endif
  146. #endif
  147.  
  148. #ifdef  CRAY
  149. #  define  Vectorize
  150. #endif
  151.  
  152. #ifdef  Vectorize
  153. #  define maxpatterns  10000  /* maximum number of different site patterns */
  154. #endif
  155.  
  156. #include <stdio.h>
  157. #include <math.h>
  158. #include "fastDNAml.h"  /*  Requires version 1.1  */
  159.  
  160. #if Master || Slave
  161. #  include "p4.h"
  162. #  include "comm_link.h" 
  163. #endif
  164.  
  165. /*  Global variables */
  166.  
  167. xarray       *usedxtip, *freextip;
  168.  
  169. #if Sequential     /*  Use standard input */
  170. #  undef   DNAML_STEP
  171. #  define  DNAML_STEP  0
  172. #ifdef NoSTDIN
  173.   FILE *INFILE;
  174. #else
  175. #  define  INFILE  stdin
  176. #endif
  177. #endif
  178.  
  179. #if Master
  180. #  define MAX_SEND_AHEAD 400
  181.   char   *best_tr_recv = NULL;     /* these are used for flow control */
  182.   double  best_lk_recv;
  183.   int     send_ahead = 0;          /* number of outstanding sends */
  184.  
  185. #  ifdef DNAML_STEP
  186. #    define  DNAML_STEP  1
  187. #  endif
  188. #  define  INFILE   Seqf
  189. #  define  OUTFILE  Outf
  190.   FILE  *INFILE, *OUTFILE;
  191.   comm_block comm_slave;
  192. #endif
  193.  
  194. #if Slave
  195. #  undef   DNAML_STEP
  196. #  define  DNAML_STEP  0
  197. #  define  INFILE   Seqf
  198. #  define  OUTFILE  Outf
  199.   FILE  *INFILE, *OUTFILE;
  200.   comm_block comm_master;
  201. #endif
  202.  
  203. #if Debug
  204.   FILE *debug;
  205. #endif
  206.  
  207. #if DNAML_STEP
  208.   int begin_step_time, end_step_time;
  209. #  define  REPORT_ADD_SPECS  p4_send(DNAML_ADD_SPECS, DNAML_HOST_ID, NULL, 0)
  210. #  define  REPORT_SEND_TREE  p4_send(DNAML_SEND_TREE, DNAML_HOST_ID, NULL, 0)
  211. #  define  REPORT_RECV_TREE  p4_send(DNAML_RECV_TREE, DNAML_HOST_ID, NULL, 0)
  212. #  define  REPORT_STEP_TIME \
  213.    {\
  214.        char send_buf[80]; \
  215.        end_step_time = p4_clock(); \
  216.        (void) sprintf(send_buf, "%d", end_step_time-begin_step_time); \
  217.        p4_send(DNAML_STEP_TIME, DNAML_HOST_ID, send_buf,strlen(send_buf)+1); \
  218.        begin_step_time = end_step_time; \
  219.    }
  220. #else
  221. #  define  REPORT_ADD_SPECS
  222. #  define  REPORT_SEND_TREE
  223. #  define  REPORT_RECV_TREE
  224. #  define  REPORT_STEP_TIME
  225. #endif
  226.  
  227. /*=======================================================================*/
  228. /*                              PROGRAM                                  */
  229. /*=======================================================================*/
  230. /*                    Best tree handling for dnaml                       */
  231. /*=======================================================================*/
  232.  
  233. /*  Tip value comparisons
  234.  *
  235.  *  Use void pointers to hide type from other routines.  Only tipValPtr and
  236.  *  cmpTipVal need to be changed to alter the nature of the values compared
  237.  *  (e.g., names instead of node numbers).
  238.  *
  239.  *    cmpTipVal(tipValPtr(nodeptr p), tipValPtr(nodeptr q)) == -1, 0 or 1.
  240.  *
  241.  *  This provides direct comparison of tip values (for example, for
  242.  *  definition of tr->start).
  243.  */
  244.  
  245.  
  246. void  *tipValPtr (p)  nodeptr  p; { return  (void *) & p->number; }
  247.  
  248.  
  249. int  cmpTipVal (v1, v2)
  250.     void  *v1, *v2;
  251.   { /* cmpTipVal */
  252.     int  i1, i2;
  253.  
  254.     i1 = *((int *) v1);
  255.     i2 = *((int *) v2);
  256.     return  (i1 < i2) ? -1 : ((i1 == i2) ? 0 : 1);
  257.   } /* cmpTipVal */
  258.  
  259.  
  260. /*  These are the only routines that need to UNDERSTAND topologies */
  261.  
  262.  
  263. topol  *setupTopol (maxtips, nsites)
  264.     int     maxtips, nsites;
  265.   { /* setupTopol */
  266.     topol   *tpl;
  267.  
  268.     if (! (tpl = (topol *) Malloc(sizeof(topol))) || 
  269.         ! (tpl->links = (connptr) Malloc((2*maxtips-3) * sizeof(connect))) || 
  270.         (nsites && ! (tpl->log_f
  271.                 = (double *) Malloc(nsites * sizeof(double))))) {
  272.       printf("ERROR: Unable to get topology memory");
  273.       tpl = (topol *) NULL;
  274.       }
  275.  
  276.     else {
  277.       if (nsites == 0)  tpl->log_f = (double *) NULL;
  278.       tpl->likelihood  = unlikely;
  279.       tpl->start       = (node *) NULL;
  280.       tpl->nextlink    = 0;
  281.       tpl->ntips       = 0;
  282.       tpl->nextnode    = 0;
  283.       tpl->opt_level   = 0;     /* degree of branch swapping explored */
  284.       tpl->scrNum      = 0;     /* position in sorted list of scores */
  285.       tpl->tplNum      = 0;     /* position in sorted list of trees */
  286.       tpl->log_f_valid = 0;     /* log_f value sites */
  287.       tpl->prelabeled  = TRUE;
  288.       tpl->smoothed    = FALSE; /* branch optimization converged? */
  289.       }
  290.  
  291.     return  tpl;
  292.   } /* setupTopol */
  293.  
  294.  
  295. void  freeTopol (tpl)
  296.     topol  *tpl;
  297.   { /* freeTopol */
  298.     Free(tpl->links);
  299.     if (tpl->log_f)  Free(tpl->log_f);
  300.     Free(tpl);
  301.   } /* freeTopol */
  302.  
  303.  
  304. int  saveSubtree (p, tpl)
  305.     nodeptr  p;
  306.     topol   *tpl;
  307.     /*  Save a subtree in a standard order so that earlier branches
  308.      *  from a node contain lower value tips than do second branches from
  309.      *  the node.  This code works with arbitrary furcations in the tree.
  310.      */
  311.   { /* saveSubtree */
  312.     connptr  r, r0;
  313.     nodeptr  q, s;
  314.     int      t, t0, t1;
  315.  
  316.     r0 = tpl->links;
  317.     r = r0 + (tpl->nextlink)++;
  318.     r->p = p;
  319.     r->q = q = p->back;
  320.     r->z = p->z;
  321.     r->descend = 0;                     /* No children (yet) */
  322.  
  323.     if (q->tip) {
  324.       r->valptr = tipValPtr(q);         /* Assign value */
  325.       }
  326.  
  327.     else {                              /* Internal node, look at children */
  328.       s = q->next;                      /* First child */
  329.       do {
  330.         t = saveSubtree(s, tpl);        /* Generate child's subtree */
  331.  
  332.         t0 = 0;                         /* Merge child into list */
  333.         t1 = r->descend;
  334.         while (t1 && (cmpTipVal(r0[t1].valptr, r0[t].valptr) < 0)) {
  335.           t0 = t1;
  336.           t1 = r0[t1].sibling;
  337.           }
  338.         if (t0) r0[t0].sibling = t;  else  r->descend = t;
  339.         r0[t].sibling = t1;
  340.  
  341.         s = s->next;                    /* Next child */
  342.         } while (s != q);
  343.  
  344.       r->valptr = r0[r->descend].valptr;   /* Inherit first child's value */
  345.       }                                 /* End of internal node processing */
  346.  
  347.     return  r - r0;
  348.   } /* saveSubtree */
  349.  
  350.  
  351. nodeptr  minSubtreeTip (p0)
  352.     nodeptr  p0;
  353.   { /* minTreeTip */
  354.     nodeptr  minTip, p, testTip;
  355.  
  356.     if (p0->tip) return p0;
  357.  
  358.     p = p0->next;
  359.     minTip = minSubtreeTip(p->back);
  360.     while ((p = p->next) != p0) {
  361.       testTip = minSubtreeTip(p->back);
  362.       if (cmpTipVal(tipValPtr(testTip), tipValPtr(minTip)) < 0)
  363.         minTip = testTip;
  364.       }
  365.     return minTip;
  366.   } /* minTreeTip */
  367.  
  368.  
  369. nodeptr  minTreeTip (p)
  370.     nodeptr  p;
  371.   { /* minTreeTip */
  372.     nodeptr  minp, minpb;
  373.  
  374.     minp  = minSubtreeTip(p);
  375.     minpb = minSubtreeTip(p->back);
  376.     return cmpTipVal(tipValPtr(minp), tipValPtr(minpb)) < 0 ? minp : minpb;
  377.   } /* minTreeTip */
  378.  
  379.  
  380. void saveTree (tr, tpl)
  381.     tree   *tr;
  382.     topol  *tpl;
  383.     /*  Save a tree topology in a standard order so that first branches
  384.      *  from a node contain lower value tips than do second branches from
  385.      *  the node.  The root tip should have the lowest value of all.
  386.      */
  387.   { /* saveTree */
  388.     connptr  r;
  389.     double  *tr_log_f, *tpl_log_f;
  390.     int  i;
  391.  
  392.     tpl->nextlink = 0;                             /* Reset link pointer */
  393.     r = tpl->links + saveSubtree(minTreeTip(tr->start), tpl);  /* Save tree */
  394.     r->sibling = 0;
  395.  
  396.     tpl->likelihood = tr->likelihood;
  397.     tpl->start      = tr->start;
  398.     tpl->ntips      = tr->ntips;
  399.     tpl->nextnode   = tr->nextnode;
  400.     tpl->opt_level  = tr->opt_level;
  401.     tpl->prelabeled = tr->prelabeled;
  402.     tpl->smoothed   = tr->smoothed;
  403.  
  404.     if (tpl_log_f = tpl->log_f) {
  405.       tr_log_f  = tr->log_f;
  406.       i = tpl->log_f_valid = tr->log_f_valid;
  407.       while (--i >= 0)  *tpl_log_f++ = *tr_log_f++;
  408.       }
  409.     else {
  410.       tpl->log_f_valid = 0;
  411.       }
  412.   } /* saveTree */
  413.  
  414.  
  415. void copyTopol (tpl1, tpl2)
  416.     topol  *tpl1, *tpl2;
  417.   { /* copyTopol */
  418.     connptr  r1, r2, r10, r20;
  419.     double  *tpl1_log_f, *tpl2_log_f;
  420.     int  i;
  421.  
  422.     r10 = tpl1->links;
  423.     r20 = tpl2->links;
  424.     tpl2->nextlink = tpl1->nextlink; 
  425.  
  426.     r1 = r10;
  427.     r2 = r20;
  428.     i = 2 * tpl1->ntips - 3;
  429.     while (--i >= 0) {
  430.       r2->z = r1->z;
  431.       r2->p = r1->p;
  432.       r2->q = r1->q;
  433.       r2->valptr = r1->valptr;
  434.       r2->descend = r1->descend; 
  435.       r2->sibling = r1->sibling; 
  436.       r1++;
  437.       r2++;
  438.       }
  439.  
  440.     if (tpl1->log_f_valid && tpl2->log_f) {
  441.       tpl1_log_f = tpl1->log_f;
  442.       tpl2_log_f = tpl2->log_f;
  443.       tpl2->log_f_valid = i = tpl1->log_f_valid;
  444.       while (--i >= 0)  *tpl2_log_f++ = *tpl1_log_f++;
  445.       }
  446.     else {
  447.       tpl2->log_f_valid = 0;
  448.       }
  449.  
  450.     tpl2->likelihood = tpl1->likelihood;
  451.     tpl2->start      = tpl1->start;
  452.     tpl2->ntips      = tpl1->ntips;
  453.     tpl2->nextnode   = tpl1->nextnode;
  454.     tpl2->opt_level  = tpl1->opt_level;
  455.     tpl2->prelabeled = tpl1->prelabeled;
  456.     tpl2->scrNum     = tpl1->scrNum;
  457.     tpl2->tplNum     = tpl1->tplNum;
  458.     tpl2->smoothed   = tpl1->smoothed;
  459.   } /* copyTopol */
  460.  
  461.  
  462. boolean restoreTree (tpl, tr)
  463.     topol  *tpl;
  464.     tree   *tr;
  465.   { /* restoreTree */
  466.     void  hookup();
  467.     boolean  initrav();
  468.  
  469.     connptr  r;
  470.     nodeptr  p, p0;
  471.     double  *tr_log_f, *tpl_log_f;
  472.     int  i;
  473.  
  474. /*  Clear existing connections */
  475.  
  476.     for (i = 1; i <= 2*(tr->mxtips) - 2; i++) {  /* Uses p = p->next at tip */
  477.       p0 = p = tr->nodep[i];
  478.       do {
  479.         p->back = (nodeptr) NULL;
  480.         p = p->next;
  481.         } while (p != p0);
  482.       }
  483.  
  484. /*  Copy connections from topology */
  485.  
  486.     for (r = tpl->links, i = 0; i < tpl->nextlink; r++, i++) {
  487.       hookup(r->p, r->q, r->z);
  488.       }
  489.  
  490.     tr->likelihood = tpl->likelihood;
  491.     tr->start      = tpl->start;
  492.     tr->ntips      = tpl->ntips;
  493.     tr->nextnode   = tpl->nextnode;
  494.     tr->opt_level  = tpl->opt_level;
  495.     tr->prelabeled = tpl->prelabeled;
  496.     tr->smoothed   = tpl->smoothed;
  497.  
  498.     if (tpl_log_f = tpl->log_f) {
  499.       tr_log_f = tr->log_f;
  500.       i = tr->log_f_valid = tpl->log_f_valid;
  501.       while (--i >= 0)  *tr_log_f++ = *tpl_log_f++;
  502.       }
  503.     else {
  504.       tr->log_f_valid = 0;
  505.       }
  506.  
  507.     return (initrav(tr, tr->start) && initrav(tr, tr->start->back));
  508.   } /* restoreTree */
  509.  
  510.  
  511. int initBestTree (bt, newkeep, numsp, sites)
  512.     bestlist  *bt;
  513.     int        newkeep, numsp, sites;
  514.   { /* initBestTree */
  515.     int  i, nlogf;
  516.  
  517.  
  518.     bt->nkeep = 0;
  519.  
  520.     if (bt->ninit <= 0) {
  521.       if (! (bt->start = setupTopol(numsp, sites)))  return  0;
  522.       bt->ninit = -1;
  523.       bt->nvalid = 0;
  524.       bt->numtrees = 0;
  525.       bt->best = unlikely;
  526.       bt->improved = FALSE;
  527.       bt->byScore = (topol **) Malloc((newkeep+1) * sizeof(topol *));
  528.       bt->byTopol = (topol **) Malloc((newkeep+1) * sizeof(topol *));
  529.       if (! bt->byScore || ! bt->byTopol) {
  530.         fprintf(stderr, "initBestTree: Malloc failure\n");
  531.         return 0;
  532.         }
  533.       }
  534.     else if (ABS(newkeep) > bt->ninit) {
  535.       if (newkeep <  0) newkeep = -(bt->ninit);
  536.       else newkeep = bt->ninit;
  537.       }
  538.  
  539.     if (newkeep < 1) {    /*  Use negative newkeep to clear list  */
  540.       newkeep = -newkeep;
  541.       if (newkeep < 1) newkeep = 1;
  542.       bt->nvalid = 0;
  543.       bt->best = unlikely;
  544.       }
  545.  
  546.     if (bt->nvalid >= newkeep) {
  547.       bt->nvalid = newkeep;
  548.       bt->worst = bt->byScore[newkeep]->likelihood;
  549.       }
  550.     else {
  551.       bt->worst = unlikely;
  552.       }
  553.  
  554.     for (i = bt->ninit + 1; i <= newkeep; i++) {
  555.       nlogf = (i <= maxlogf) ? sites : 0;
  556.       if (! (bt->byScore[i] = setupTopol(numsp, nlogf)))  break;
  557.       bt->byTopol[i] = bt->byScore[i];
  558.       bt->ninit = i;
  559.       }
  560.  
  561.     return  (bt->nkeep = MIN(newkeep, bt->ninit));
  562.   } /* initBestTree */
  563.  
  564.  
  565.  
  566. int resetBestTree (bt)
  567.     bestlist  *bt;
  568.   { /* resetBestTree */
  569.     bt->best = unlikely;
  570.     bt->worst = unlikely;
  571.     bt->nvalid = 0;
  572.     bt->improved = FALSE;
  573.   } /* resetBestTree */
  574.  
  575.  
  576. boolean  freeBestTree(bt)
  577.     bestlist  *bt;
  578.   { /* freeBestTree */
  579.     while (bt->ninit >= 0)  freeTopol(bt->byScore[(bt->ninit)--]);
  580.     freeTopol(bt->start);
  581.     return TRUE;
  582.   } /* freeBestTree */
  583.  
  584.  
  585. /*  Compare two trees, assuming that each is in standard order.  Return
  586.  *  -1 if first preceeds second, 0 if they are identical, or +1 if first
  587.  *  follows second in standard order.  Lower number tips preceed higher
  588.  *  number tips.  A tip preceeds a corresponding internal node.  Internal
  589.  *  nodes are ranked by their lowest number tip.
  590.  */
  591.  
  592. int  cmpSubtopol (p10, p1, p20, p2)
  593.     connptr  p10, p1, p20, p2;
  594.   { /* cmpSubtopol */
  595.     connptr  p1d, p2d;
  596.     int  cmp;
  597.  
  598.     if (! p1->descend && ! p2->descend)          /* Two tips */
  599.       return cmpTipVal(p1->valptr, p2->valptr);
  600.  
  601.     if (! p1->descend) return -1;                /* p1 = tip, p2 = node */
  602.     if (! p2->descend) return  1;                /* p2 = tip, p1 = node */
  603.  
  604.     p1d = p10 + p1->descend;
  605.     p2d = p20 + p2->descend;
  606.     while (1) {                                  /* Two nodes */
  607.       if (cmp = cmpSubtopol(p10, p1d, p20, p2d))  return cmp; /* Subtrees */
  608.       if (! p1d->sibling && ! p2d->sibling)  return  0; /* Lists done */
  609.       if (! p1d->sibling) return -1;             /* One done, other not */
  610.       if (! p2d->sibling) return  1;             /* One done, other not */
  611.       p1d = p10 + p1d->sibling;                  /* Neither done */
  612.       p2d = p20 + p2d->sibling;
  613.       }
  614.   } /* cmpSubtopol */
  615.  
  616.  
  617.  
  618. int  cmpTopol (tpl1, tpl2)
  619.     void  *tpl1, *tpl2;
  620.   { /* cmpTopol */
  621.     connptr  r1, r2;
  622.     int      cmp;
  623.  
  624.     r1 = ((topol *) tpl1)->links;
  625.     r2 = ((topol *) tpl2)->links;
  626.     cmp = cmpTipVal(tipValPtr(r1->p), tipValPtr(r2->p));
  627.     if (cmp) return cmp;
  628.     return  cmpSubtopol(r1, r1, r2, r2);
  629.   } /* cmpTopol */
  630.  
  631.  
  632.  
  633. int  cmpTplScore (tpl1, tpl2)
  634.     void  *tpl1, *tpl2;
  635.   { /* cmpTplScore */
  636.     double  l1, l2;
  637.  
  638.     l1 = ((topol *) tpl1)->likelihood;
  639.     l2 = ((topol *) tpl2)->likelihood;
  640.     return  (l1 > l2) ? -1 : ((l1 == l2) ? 0 : 1);
  641.   } /* cmpTplScore */
  642.  
  643.  
  644.  
  645. /*  Find an item in a sorted list of n items.  If the item is in the list,
  646.  *  return its index.  If it is not in the list, return the negative of the
  647.  *  position into which it should be inserted.
  648.  */
  649.  
  650. int  findInList (item, list, n, cmpFunc)
  651.     void  *item;
  652.     void  *list[];
  653.     int    n;
  654.     int  (* cmpFunc)();
  655.   { /* findInList */
  656.     int  mid, hi, lo, cmp;
  657.  
  658.     if (n < 1) return  -1;                    /*  No match; first index  */
  659.  
  660.     lo = 1;
  661.     mid = 0;
  662.     hi = n;
  663.     while (lo < hi) {
  664.       mid = (lo + hi) >> 1;
  665.       cmp = (* cmpFunc)(item, list[mid-1]);
  666.       if (cmp) {
  667.         if (cmp < 0) hi = mid;
  668.         else lo = mid + 1;
  669.         }
  670.       else  return  mid;                        /*  Exact match  */
  671.       }
  672.  
  673.     if (lo != mid) {
  674.        cmp = (* cmpFunc)(item, list[lo-1]);
  675.        if (cmp == 0) return lo;
  676.        }
  677.     if (cmp > 0) lo++;                         /*  Result of step = 0 test  */
  678.     return  -lo;
  679.   } /* findInList */
  680.  
  681.  
  682.  
  683. int  findTreeInList (bt, tr)
  684.     bestlist  *bt;
  685.     tree      *tr;
  686.   { /* findTreeInList */
  687.     topol  *tpl;
  688.  
  689.     tpl = bt->byScore[0];
  690.     saveTree(tr, tpl);
  691.     return  findInList((void *) tpl, (void **) (& (bt->byTopol[1])),
  692.                        bt->nvalid, cmpTopol);
  693.   } /* findTreeInList */
  694.  
  695.  
  696. int  saveBestTree (bt, tr)
  697.     bestlist  *bt;
  698.     tree      *tr;
  699.   { /* saveBestTree */
  700.     double *tr_log_f, *tpl_log_f;
  701.     topol  *tpl, *reuse;
  702.     int  tplNum, scrNum, reuseScrNum, reuseTplNum, i, oldValid, newValid;
  703.  
  704.     tplNum = findTreeInList(bt, tr);
  705.     tpl = bt->byScore[0];
  706.     oldValid = newValid = bt->nvalid;
  707.  
  708.     if (tplNum > 0) {                      /* Topology is in list  */
  709.       reuse = bt->byTopol[tplNum];         /* Matching topol  */
  710.       reuseScrNum = reuse->scrNum;
  711.       reuseTplNum = reuse->tplNum;
  712.       }
  713.                                            /* Good enough to keep? */
  714.     else if (tr->likelihood < bt->worst)  return 0;
  715.  
  716.     else {                                 /* Topology is not in list */
  717.       tplNum = -tplNum;                    /* Add to list (not replace) */
  718.       if (newValid < bt->nkeep) bt->nvalid = ++newValid;
  719.       reuseScrNum = newValid;              /* Take worst tree */
  720.       reuse = bt->byScore[reuseScrNum];
  721.       reuseTplNum = (newValid > oldValid) ? newValid : reuse->tplNum;
  722.       if (tr->likelihood > bt->start->likelihood) bt->improved = TRUE;
  723.       }
  724.  
  725.     scrNum = findInList((void *) tpl, (void **) (& (bt->byScore[1])),
  726.                          oldValid, cmpTplScore);
  727.     scrNum = ABS(scrNum);
  728.  
  729.     if (scrNum < reuseScrNum)
  730.       for (i = reuseScrNum; i > scrNum; i--)
  731.         (bt->byScore[i] = bt->byScore[i-1])->scrNum = i;
  732.  
  733.     else if (scrNum > reuseScrNum) {
  734.       scrNum--;
  735.       for (i = reuseScrNum; i < scrNum; i++)
  736.         (bt->byScore[i] = bt->byScore[i+1])->scrNum = i;
  737.       }
  738.  
  739.     if (tplNum < reuseTplNum)
  740.       for (i = reuseTplNum; i > tplNum; i--)
  741.         (bt->byTopol[i] = bt->byTopol[i-1])->tplNum = i;
  742.  
  743.     else if (tplNum > reuseTplNum) {
  744.       tplNum--;
  745.       for (i = reuseTplNum; i < tplNum; i++)
  746.         (bt->byTopol[i] = bt->byTopol[i+1])->tplNum = i;
  747.       }
  748.  
  749.     if (tpl_log_f = tpl->log_f) {
  750.       tr_log_f  = tr->log_f;
  751.       i = tpl->log_f_valid = tr->log_f_valid;
  752.       while (--i >= 0)  *tpl_log_f++ = *tr_log_f++;
  753.       }
  754.     else {
  755.       tpl->log_f_valid = 0;
  756.       }
  757.  
  758.     tpl->scrNum = scrNum;
  759.     tpl->tplNum = tplNum;
  760.     bt->byTopol[tplNum] = bt->byScore[scrNum] = tpl;
  761.     bt->byScore[0] = reuse;
  762.  
  763.     if (scrNum == 1)  bt->best = tr->likelihood;
  764.     if (newValid == bt->nkeep) bt->worst = bt->byScore[newValid]->likelihood;
  765.  
  766.     return  scrNum;
  767.   } /* saveBestTree */
  768.  
  769.  
  770. int  startOpt (bt, tr)
  771.     bestlist  *bt;
  772.     tree      *tr;
  773.   { /* startOpt */
  774.     int  scrNum;
  775.  
  776.     scrNum = saveBestTree(bt, tr);
  777.     copyTopol(bt->byScore[scrNum], bt->start);
  778.     bt->improved = FALSE;
  779.     return  scrNum;
  780.   } /* startOpt */
  781.  
  782.  
  783. int  setOptLevel (bt, opt_level)
  784.     bestlist *bt;
  785.     int       opt_level;
  786.   { /* setOptLevel */
  787.     int  tplNum, scrNum;
  788.  
  789.     tplNum = findInList((void *) bt->start, (void **) (&(bt->byTopol[1])),
  790.                         bt->nvalid, cmpTopol);
  791.     if (tplNum > 0) {
  792.       bt->byTopol[tplNum]->opt_level = opt_level;
  793.       scrNum = bt->byTopol[tplNum]->scrNum;
  794.       }
  795.     else {
  796.       scrNum = 0;
  797.       }
  798.  
  799.     return  scrNum;
  800.   } /* setOptLevel */
  801.  
  802.  
  803. int  recallBestTree (bt, rank, tr)
  804.     bestlist  *bt;
  805.     int        rank;
  806.     tree      *tr;
  807.   { /* recallBestTree */
  808.     if (rank < 1)  rank = 1;
  809.     if (rank > bt->nvalid)  rank = bt->nvalid;
  810.     if (rank > 0)  if (! restoreTree(bt->byScore[rank], tr)) return FALSE;
  811.     return  rank;
  812.   } /* recallBestTree */
  813.  
  814.  
  815. /*=======================================================================*/
  816. /*                       End of best tree routines                       */
  817. /*=======================================================================*/
  818.  
  819.  
  820. #if 0
  821.   void hang(msg) char *msg; {printf("Hanging around: %s\n", msg); while(1);}
  822. #endif
  823.  
  824.  
  825. boolean getnums (rdta)
  826.     rawDNA   *rdta;
  827.     /* input number of species, number of sites */
  828.   { /* getnums */
  829.     printf("\n%s, version %s, %s\n\n",
  830.             programName,
  831.             programVersion,
  832.             programDate);
  833.     printf("Based on Joseph Felsenstein's\n\n");
  834.     printf("   Nucleic acid sequence Maximum Likelihood method, ");
  835.     printf("version 3.3\n\n");
  836.  
  837.     if (fscanf(INFILE, "%d %d", & rdta->numsp, & rdta->sites) != 2) {
  838.       printf("ERROR: Problem reading number of species and sites\n");
  839.       return FALSE;
  840.       }
  841.     printf("%d Species, %d Sites\n\n", rdta->numsp, rdta->sites);
  842.  
  843.     if (rdta->numsp < 4) {
  844.       printf("TOO FEW SPECIES\n");
  845.       return FALSE;
  846.       }
  847.  
  848.     if (rdta->sites < 1) {
  849.       printf("TOO FEW SITES\n");
  850.       return FALSE;
  851.       }
  852.  
  853.     return TRUE;
  854.   } /* getnums */
  855.  
  856.  
  857. boolean digitchar (ch) int ch; {return (ch >= '0' && ch <= '9'); }
  858.  
  859.  
  860. boolean whitechar (ch) int ch;
  861.    { return (ch == ' ' || ch == '\n' || ch == '\t');
  862.      }
  863.  
  864.  
  865. void uppercase (chptr)
  866.       int *chptr;
  867.     /* convert character to upper case -- either ASCII or EBCDIC */
  868.   { /* uppercase */
  869.     int  ch;
  870.  
  871.     ch = *chptr;
  872.     if ((ch >= 'a' && ch <= 'i') || (ch >= 'j' && ch <= 'r')
  873.                                  || (ch >= 's' && ch <= 'z'))
  874.       *chptr = ch + 'A' - 'a';
  875.   } /* uppercase */
  876.  
  877.  
  878. int base36 (ch)
  879.     int ch;
  880.   { /* base36 */
  881.     if      (ch >= '0' && ch <= '9') return (ch - '0');
  882.     else if (ch >= 'A' && ch <= 'I') return (ch - 'A' + 10);
  883.     else if (ch >= 'J' && ch <= 'R') return (ch - 'J' + 19);
  884.     else if (ch >= 'S' && ch <= 'Z') return (ch - 'S' + 28);
  885.     else if (ch >= 'a' && ch <= 'i') return (ch - 'a' + 10);
  886.     else if (ch >= 'j' && ch <= 'r') return (ch - 'j' + 19);
  887.     else if (ch >= 's' && ch <= 'z') return (ch - 's' + 28);
  888.     else return -1;
  889.   } /* base36 */
  890.  
  891.  
  892. int itobase36 (i)
  893.     int  i;
  894.   { /* itobase36 */
  895.     if      (i <  0) return '?';
  896.     else if (i < 10) return (i      + '0');
  897.     else if (i < 19) return (i - 10 + 'A');
  898.     else if (i < 28) return (i - 19 + 'J');
  899.     else if (i < 36) return (i - 28 + 'S');
  900.     else return '?';
  901.   } /* itobase36 */
  902.  
  903.  
  904. int findch (c)
  905.     int  c;
  906.   { /* findch */
  907.     int ch;
  908.  
  909.     while ((ch = getc(INFILE)) != EOF && ch != c) ;
  910.     return  ch;
  911.   } /* findch */
  912.  
  913.  
  914. #if Master || Slave
  915. int str_findch (treestrp, c)
  916.     char **treestrp;
  917.     int  c;
  918.   { /* str_findch */
  919.     int ch;
  920.  
  921.     while ((ch = *(*treestrp)++) != NULL && ch != c) ;
  922.     return  ch;
  923.   } /* str_findch */
  924. #endif
  925.  
  926.  
  927. boolean inputboot(adef)
  928.     analdef      *adef;
  929.     /* read the bootstrap auxilliary info */
  930.   { /* inputboot */
  931.     if (! adef->boot) {
  932.       printf("ERROR: Unexpected Bootstrap auxiliary data line\n");
  933.       return FALSE;
  934.       }
  935.     else if (fscanf(INFILE, "%ld", & adef->boot) != 1 ||
  936.              findch('\n') == EOF) {
  937.       printf("ERROR: Problem reading boostrap random seed value\n");
  938.       return FALSE;
  939.       }
  940.  
  941.     return TRUE;
  942.   } /* inputboot */
  943.  
  944.  
  945. boolean inputcategories (rdta)
  946.     rawDNA       *rdta;
  947.     /* read the category rates and the categories for each site */
  948.   { /* inputcategories */
  949.     int  i, j, ch, ci;
  950.  
  951.     if (rdta->categs >= 0) {
  952.       printf("ERROR: Unexpected Categories auxiliary data line\n");
  953.       return FALSE;
  954.       }
  955.     if (fscanf(INFILE, "%d", & rdta->categs) != 1) {
  956.       printf("ERROR: Problem reading number of rate categories\n");
  957.       return FALSE;
  958.       }
  959.     if (rdta->categs < 1 || rdta->categs > maxcategories) {
  960.       printf("ERROR: Bad number of categories: %d\n", rdta->categs);
  961.       printf("Must be in range 1 - %d\n", maxcategories);
  962.       return FALSE;
  963.       }
  964.  
  965.     for (j = 1; j <= rdta->categs
  966.            && fscanf(INFILE, "%lf", &(rdta->catrat[j])) == 1; j++) ;
  967.  
  968.     if ((j <= rdta->categs) || (findch('\n') == EOF)) {
  969.       printf("ERROR: Problem reading rate values\n");
  970.       return FALSE;
  971.       }
  972.  
  973.     for (i = 1; i <= nmlngth; i++)  (void) getc(INFILE);
  974.  
  975.     i = 1;
  976.     while (i <= rdta->sites) {
  977.       ch = getc(INFILE);
  978.       ci = base36(ch);
  979.       if (ci >= 0 && ci <= rdta->categs)
  980.         rdta->sitecat[i++] = ci;
  981.       else if (! whitechar(ch)) {
  982.         printf("ERROR: Bad category character (%c) at site %d\n", ch, i);
  983.         return FALSE;
  984.         }
  985.       }
  986.  
  987.     if (findch('\n') == EOF) {      /* skip to end of line */
  988.       printf("ERROR: Missing newline at end of category data\n");
  989.       return FALSE;
  990.       }
  991.  
  992.     return TRUE;
  993.   } /* inputcategories */
  994.  
  995.  
  996. boolean inputextra (adef)
  997.     analdef  *adef;
  998.   { /* inputextra */
  999.     if (fscanf(INFILE,"%d", & adef->extra) != 1 ||
  1000.         findch('\n') == EOF) {
  1001.       printf("ERROR: Problem reading extra info value\n");
  1002.       return FALSE;
  1003.       }
  1004.  
  1005.     return TRUE;
  1006.   } /* inputextra */
  1007.  
  1008.  
  1009. boolean inputfreqs (rdta)
  1010.     rawDNA   *rdta;
  1011.   { /* inputfreqs */
  1012.     if (fscanf(INFILE, "%lf%lf%lf%lf",
  1013.                & rdta->freqa, & rdta->freqc,
  1014.                & rdta->freqg, & rdta->freqt) != 4 ||
  1015.         findch('\n') == EOF) {
  1016.       printf("ERROR: Problem reading user base frequencies data\n");
  1017.       return  FALSE;
  1018.       }
  1019.  
  1020.     rdta->freqread = TRUE;
  1021.     return TRUE;
  1022.   } /* inputfreqs */
  1023.  
  1024.  
  1025. boolean inputglobal (tr)
  1026.     tree     *tr;
  1027.     /* input the global option information */
  1028.   { /* inputglobal */
  1029.     int  ch;
  1030.  
  1031.     if (tr->global != -2) {
  1032.       printf("ERROR: Unexpected Global auxiliary data line\n");
  1033.       return FALSE;
  1034.       }
  1035.     if (fscanf(INFILE, "%d", &(tr->global)) != 1) {
  1036.       printf("ERROR: Problem reading rearrangement region size\n");
  1037.       return FALSE;
  1038.       }
  1039.     if (tr->global < 0) {
  1040.       printf("WARNING: Global region size too small;\n");
  1041.       printf("         value reset to local\n\n");
  1042.       tr->global = 1;
  1043.       }
  1044.     else if (tr->global == 0)  tr->partswap = 0;
  1045.     else if (tr->global > tr->mxtips - 3) {
  1046.       tr->global = tr->mxtips - 3;
  1047.       }
  1048.  
  1049.     while ((ch = getc(INFILE)) != '\n') {  /* Scan for second value */
  1050.       if (! whitechar(ch)) {
  1051.         if (ch != EOF)  (void) ungetc(ch, INFILE);
  1052.         if (ch == EOF || fscanf(INFILE, "%d", &(tr->partswap)) != 1
  1053.                       || findch('\n') == EOF) {
  1054.           printf("ERROR: Problem reading insert swap region size\n");
  1055.           return FALSE;
  1056.           }
  1057.         else if (tr->partswap < 0)  tr->partswap = 1;
  1058.         else if (tr->partswap > tr->mxtips - 3) {
  1059.           tr->partswap = tr->mxtips - 3;
  1060.           }
  1061.  
  1062.         if (tr->partswap > tr->global)  tr->global = tr->partswap;
  1063.         break;   /*  Break while loop */
  1064.         }
  1065.       }
  1066.  
  1067.     return TRUE;
  1068.   } /* inputglobal */
  1069.  
  1070.  
  1071. boolean inputjumble (adef)
  1072.     analdef  *adef;
  1073.   { /* inputjumble */
  1074.     if (! adef->jumble) {
  1075.       printf("ERROR: Unexpected Jumble auxiliary data line\n");
  1076.       return FALSE;
  1077.       }
  1078.     else if (fscanf(INFILE, "%ld", & adef->jumble) != 1 ||
  1079.              findch('\n') == EOF) {
  1080.       printf("ERROR: Problem reading jumble random seed value\n");
  1081.       return FALSE;
  1082.       }
  1083.     else if (adef->jumble == 0) {
  1084.       printf("WARNING: Jumble random number seed is zero\n\n");
  1085.       }
  1086.  
  1087.     return TRUE;
  1088.   } /* inputjumble */
  1089.  
  1090.  
  1091. boolean inputkeep (adef)
  1092.     analdef  *adef;
  1093.   { /* inputkeep */
  1094.     if (fscanf(INFILE, "%d", & adef->nkeep) != 1 ||
  1095.         findch('\n') == EOF || adef->nkeep < 1) {
  1096.       printf("ERROR: Problem reading number of kept trees\n");
  1097.       return FALSE;
  1098.       }
  1099.  
  1100.     return TRUE;
  1101.   } /* inputkeep */
  1102.  
  1103.  
  1104. boolean inputoutgroup (adef, tr)
  1105.     analdef  *adef;
  1106.     tree     *tr;
  1107.   { /* inputoutgroup */
  1108.     if (! adef->root || tr->outgr > 0) {
  1109.       printf("ERROR: Unexpected Outgroup auxiliary data line\n");
  1110.       return FALSE;
  1111.       }
  1112.     else if (fscanf(INFILE, "%d", &(tr->outgr)) != 1 ||
  1113.              findch('\n') == EOF) {
  1114.       printf("ERROR: Problem reading outgroup number\n");
  1115.       return FALSE;
  1116.       }
  1117.     else if ((tr->outgr < 1) || (tr->outgr > tr->mxtips)) {
  1118.       printf("ERROR: Bad outgroup: '%d'\n", tr->outgr);
  1119.       return FALSE;
  1120.       }
  1121.  
  1122.     return TRUE;
  1123.   } /* inputoutgroup */
  1124.  
  1125.  
  1126. boolean inputratio (rdta)
  1127.     rawDNA   *rdta;
  1128.   { /* inputratio */
  1129.     if (rdta->ttratio >= 0.0) {
  1130.       printf("ERROR: Unexpected Transition/transversion auxiliary data\n");
  1131.       return FALSE;
  1132.       }
  1133.     else if (fscanf(INFILE,"%lf", & rdta->ttratio)!=1 ||
  1134.              findch('\n') == EOF) {
  1135.       printf("ERROR: Problem reading transition/transversion ratio\n");
  1136.       return FALSE;
  1137.       }
  1138.  
  1139.     return TRUE;
  1140.   } /* inputratio */
  1141.  
  1142.  
  1143. /*  Y 0 is treeNone   (no tree)
  1144.     Y 1 is treeNewick
  1145.     Y 2 is treeProlog
  1146.     Y 3 is treePHYLIP
  1147.  */
  1148.  
  1149. boolean inputtreeopt (adef)
  1150.     analdef  *adef;
  1151.   { /* inputtreeopt */
  1152.     if (! adef->trout) {
  1153.       printf("ERROR: Unexpected Treefile auxiliary data\n");
  1154.       return FALSE;
  1155.       }
  1156.     else if (fscanf(INFILE,"%d", & adef->trout) != 1 ||
  1157.              findch('\n') == EOF) {
  1158.       printf("ERROR: Problem reading output tree-type number\n");
  1159.       return FALSE;
  1160.       }
  1161.     else if ((adef->trout < 0) || (adef->trout > treeMaxType)) {
  1162.       printf("ERROR: Bad output tree-type number: '%d'\n", adef->trout);
  1163.       return FALSE;
  1164.       }
  1165.  
  1166.     return TRUE;
  1167.   } /* inputtreeopt */
  1168.  
  1169.  
  1170. boolean inputweights (adef, rdta, cdta)
  1171.     analdef      *adef;
  1172.     rawDNA       *rdta;
  1173.     crunchedDNA  *cdta;
  1174.     /* input the character weights 0, 1, 2 ... 9, A, B, ... Y, Z */
  1175.   { /* inputweights */
  1176.     int i, ch, wi;
  1177.  
  1178.     if (! adef->userwgt || cdta->wgtsum > 0) {
  1179.       printf("ERROR: Unexpected Weights auxiliary data\n");
  1180.       return FALSE;
  1181.       }
  1182.  
  1183.     for (i = 2; i <= nmlngth; i++)  (void) getc(INFILE);
  1184.     cdta->wgtsum = 0;
  1185.     i = 1;
  1186.     while (i <= rdta->sites) {
  1187.       ch = getc(INFILE);
  1188.       wi = base36(ch);
  1189.       if (wi >= 0)
  1190.         cdta->wgtsum += rdta->wgt[i++] = wi;
  1191.       else if (! whitechar(ch)) {
  1192.         printf("ERROR: Bad weight character: '%c'", ch);
  1193.         printf("       Weights in dnaml must be a digit or a letter.\n");
  1194.         return FALSE;
  1195.         }
  1196.       }
  1197.  
  1198.     if (findch('\n') == EOF) {      /* skip to end of line */
  1199.       printf("ERROR: Missing newline at end of weight data\n");
  1200.       return FALSE;
  1201.       }
  1202.  
  1203.     return TRUE;
  1204.   } /* inputweights */
  1205.  
  1206.  
  1207. boolean getoptions (adef, rdta, cdta, tr)
  1208.     analdef      *adef;
  1209.     rawDNA       *rdta;
  1210.     crunchedDNA  *cdta;
  1211.     tree         *tr;
  1212.   { /* getoptions */
  1213.     int     ch, i, extranum;
  1214.  
  1215.     adef->boot    =           0;  /* Don't bootstrap column weights */
  1216.     adef->empf    =        TRUE;  /* Use empirical base frequencies */
  1217.     adef->extra   =           0;  /* No extra runtime info unless requested */
  1218.     adef->interleaved =    TRUE;  /* By default, data format is interleaved */
  1219.     adef->jumble  =       FALSE;  /* Use random addition sequence */
  1220.     adef->nkeep   =           0;  /* Keep only the one best tree */
  1221.     adef->prdata  =       FALSE;  /* Don't echo data to output stream */
  1222.     adef->qadd    =        TRUE;  /* Smooth branches globally in add */
  1223.     adef->restart =       FALSE;  /* Restart from user tree */
  1224.     adef->root    =       FALSE;  /* User-defined outgroup rooting */
  1225.     adef->trout   = treeDefType;  /* Output tree file */
  1226.     adef->trprint =        TRUE;  /* Print tree to output stream */
  1227.     rdta->categs  =           0;  /* No rate categories */
  1228.     rdta->catrat[1] =       1.0;  /* Rate values */
  1229.     rdta->freqread =      FALSE;  /* User-defined frequencies not read yet */
  1230.     rdta->ttratio =         2.0;  /* Transition/transversion rate ratio */
  1231.     tr->global    =          -1;  /* Default search locale for optimum */
  1232.     tr->mxtips    = rdta->numsp;
  1233.     tr->outgr     =           1;  /* Outgroup number */
  1234.     tr->partswap  =           1;  /* Default to swap locally after insert */
  1235.     tr->userlen   =       FALSE;  /* User-supplied branch lengths */
  1236.     adef->usertree =      FALSE;  /* User-defined tree topologies */
  1237.     adef->userwgt =       FALSE;  /* User-defined position weights */
  1238.     extranum      =           0;
  1239.  
  1240.     while ((ch = getc(INFILE)) != '\n' && ch != EOF) {
  1241.       uppercase(& ch);
  1242.       switch (ch) {
  1243.           case '1' : adef->prdata  = ! adef->prdata; break;
  1244.           case '3' : adef->trprint = ! adef->trprint; break;
  1245.           case '4' : adef->trout   = treeDefType - adef->trout; break;
  1246.           case 'B' : adef->boot    =  1; extranum++; break;
  1247.           case 'C' : rdta->categs  = -1; extranum++; break;
  1248.           case 'E' : adef->extra   = -1; break;
  1249.           case 'F' : adef->empf    = ! adef->empf; break;
  1250.           case 'G' : tr->global    = -2; break;
  1251.           case 'I' : adef->interleaved = ! adef->interleaved; break;
  1252.           case 'J' : adef->jumble  = 1; extranum++; break;
  1253.           case 'K' : extranum++; break;
  1254.           case 'L' : tr->userlen   = TRUE; break;
  1255.           case 'O' : adef->root    = TRUE; tr->outgr = 0; extranum++; break;
  1256.           case 'Q' : adef->qadd    = FALSE; break;
  1257.           case 'R' : adef->restart = TRUE; break;
  1258.           case 'T' : rdta->ttratio = -1.0; extranum++; break;
  1259.           case 'U' : adef->usertree = TRUE; break;
  1260.           case 'W' : adef->userwgt = TRUE; cdta->wgtsum = 0; extranum++; break;
  1261.           case 'Y' : adef->trout   = treeDefType - adef->trout; break;
  1262.           case ' ' : break;
  1263.           case '\t': break;
  1264.           default  :
  1265.               printf("ERROR: Bad option character: '%c'\n", ch);
  1266.               return FALSE;
  1267.           }
  1268.       }
  1269.  
  1270.     if (ch == EOF) {
  1271.       printf("ERROR: End-of-file in options list\n");
  1272.       return FALSE;
  1273.       }
  1274.  
  1275.     if (adef->usertree && adef->restart) {
  1276.       printf("ERROR:  The restart and user-tree options conflict:\n");
  1277.       printf("        Restart adds rest of taxa to a starting tree;\n");
  1278.       printf("        User-tree does not add any taxa.\n\n");
  1279.       return FALSE;
  1280.       }
  1281.  
  1282.     if (adef->usertree && adef->jumble) {
  1283.       printf("WARNING:  The jumble and user-tree options conflict:\n");
  1284.       printf("          Jumble adds taxa to a tree in random order;\n");
  1285.       printf("          User-tree does not use taxa addition.\n");
  1286.       printf("          Jumble option cancelled for this run.\n\n");
  1287.       adef->jumble = FALSE;
  1288.       }
  1289.  
  1290.     if (tr->userlen && tr->global != -1) {
  1291.       printf("ERROR:  The global and user-lengths options conflict:\n");
  1292.       printf("        Global optimizes a starting tree;\n");
  1293.       printf("        User-lengths constrain the starting tree.\n\n");
  1294.       return FALSE;
  1295.       }
  1296.  
  1297.     if (tr->userlen && ! adef->usertree) {
  1298.       printf("WARNING:  User lengths required user tree option.\n");
  1299.       printf("          User-tree option set for this run.\n\n");
  1300.       adef->usertree = TRUE;
  1301.       }
  1302.  
  1303.     rdta->wgt      = (int *)    Malloc((rdta->sites + 1) * sizeof(int));
  1304.     rdta->wgt2     = (int *)    Malloc((rdta->sites + 1) * sizeof(int));
  1305.     rdta->sitecat  = (int *)    Malloc((rdta->sites + 1) * sizeof(int));
  1306.     cdta->alias    = (int *)    Malloc((rdta->sites + 1) * sizeof(int));
  1307.     cdta->aliaswgt = (int *)    Malloc((rdta->sites + 1) * sizeof(int));
  1308.     cdta->patcat   = (int *)    Malloc((rdta->sites + 1) * sizeof(int));
  1309.     cdta->patrat   = (double *) Malloc((rdta->sites + 1) * sizeof(double));
  1310.     cdta->wr       = (double *) Malloc((rdta->sites + 1) * sizeof(double));
  1311.     cdta->wr2      = (double *) Malloc((rdta->sites + 1) * sizeof(double));
  1312.     if ( ! rdta->wgt || ! rdta->wgt2     || ! rdta->sitecat || ! cdta->alias
  1313.                      || ! cdta->aliaswgt || ! cdta->patcat  || ! cdta->patrat
  1314.                      || ! cdta->wr       || ! cdta->wr2) {
  1315.       fprintf(stderr, "getoptions: Malloc failure\n");
  1316.       return 0;
  1317.       }
  1318.  
  1319.     /*  process lines with auxiliary data */
  1320.  
  1321.     while (extranum--) {
  1322.       ch = getc(INFILE);
  1323.       uppercase(& ch);
  1324.       switch (ch) {
  1325.         case 'B':  if (! inputboot(adef)) return FALSE; break;
  1326.         case 'C':  if (! inputcategories(rdta)) return FALSE; break;
  1327.         case 'E':  if (! inputextra(adef)) return FALSE; extranum++; break;
  1328.         case 'F':  if (! inputfreqs(rdta)) return FALSE; break;
  1329.         case 'G':  if (! inputglobal(tr)) return FALSE; extranum++; break;
  1330.         case 'J':  if (! inputjumble(adef)) return FALSE; break;
  1331.         case 'K':  if (! inputkeep(adef)) return FALSE; break;
  1332.         case 'O':  if (! inputoutgroup(adef, tr)) return FALSE; break;
  1333.         case 'T':  if (! inputratio(rdta)) return FALSE; break;
  1334.         case 'W':  if (! inputweights(adef, rdta, cdta)) return FALSE; break;
  1335.         case 'Y':  if (! inputtreeopt(adef)) return FALSE; extranum++; break;
  1336.  
  1337.         default:
  1338.             printf("ERROR: Auxiliary options line starts with '%c'\n", ch);
  1339.             return FALSE;
  1340.         }
  1341.       }
  1342.  
  1343.     if (! adef->userwgt) {
  1344.       for (i = 1; i <= rdta->sites; i++) rdta->wgt[i] = 1;
  1345.       cdta->wgtsum = rdta->sites;
  1346.       }
  1347.  
  1348.     if (adef->userwgt && cdta->wgtsum < 1) {
  1349.       printf("ERROR:  Missing or bad user-supplied weight data.\n");
  1350.       return FALSE;
  1351.       }
  1352.  
  1353.     if (adef->boot) {
  1354.       printf("Bootstrap random number seed = %ld\n\n", adef->boot);
  1355.       }
  1356.  
  1357.     if (adef->jumble) {
  1358.       printf("Jumble random number seed = %ld\n\n", adef->jumble);
  1359.       }
  1360.  
  1361.     if (adef->qadd) {
  1362.       printf("Quick add (only local branches initially optimized) in effect\n\n");
  1363.       }
  1364.  
  1365.     if (rdta->categs > 0) {
  1366.       printf("Site category   Rate of change\n\n");
  1367.       for (i = 1; i <= rdta->categs; i++)
  1368.         printf("           %c%13.3f\n", itobase36(i), rdta->catrat[i]);
  1369.       putchar('\n');
  1370.       for (i = 1; i <= rdta->sites; i++) {
  1371.         if ((rdta->wgt[i] > 0) && (rdta->sitecat[i] < 1)) {
  1372.           printf("ERROR: Bad category (%c) at site %d\n",
  1373.                   itobase36(rdta->sitecat[i]), i);
  1374.           return FALSE;
  1375.           }
  1376.         }
  1377.       }
  1378.     else if (rdta->categs < 0) {
  1379.       printf("ERROR: Category auxiliary data missing from input\n");
  1380.       return FALSE;
  1381.       }
  1382.     else {                                        /* rdta->categs == 0 */
  1383.       for (i = 1; i <= rdta->sites; i++) rdta->sitecat[i] = 1;
  1384.       rdta->categs = 1;
  1385.       }
  1386.  
  1387.     if (tr->outgr < 1) {
  1388.       printf("ERROR: Outgroup auxiliary data missing from input\n");
  1389.       return FALSE;
  1390.       }
  1391.  
  1392.     if (rdta->ttratio < 0.0) {
  1393.       printf("ERROR: Transition/transversion auxiliary data missing from input\n");
  1394.       return FALSE;
  1395.       }
  1396.  
  1397.     if (tr->global < 0) {
  1398.       if (tr->global == -2) tr->global = tr->mxtips - 3;  /* Default global */
  1399.       else                  tr->global = adef->usertree ? 0 : 1;/* No global */
  1400.       }
  1401.  
  1402.     if (adef->restart) {
  1403.       printf("Restart option in effect.  ");
  1404.       printf("Sequence addition will start from appended tree.\n\n");
  1405.       }
  1406.  
  1407.     if (adef->usertree && ! tr->global) {
  1408.       printf("User-supplied tree topology%swill be used.\n\n",
  1409.         tr->userlen ? " and branch lengths " : " ");
  1410.       }
  1411.     else {
  1412.       if (! adef->usertree) {
  1413.         printf("Rearrangements of partial trees may cross %d %s.\n",
  1414.                tr->partswap, tr->partswap == 1 ? "branch" : "branches");
  1415.         }
  1416.       printf("Rearrangements of full tree may cross %d %s.\n\n",
  1417.              tr->global, tr->global == 1 ? "branch" : "branches");
  1418.       }
  1419.  
  1420.     if (! adef->usertree && adef->nkeep == 0) adef->nkeep = 1;
  1421.  
  1422.     return TRUE;
  1423.   } /* getoptions */
  1424.  
  1425.  
  1426. boolean getbasefreqs (rdta)
  1427.     rawDNA   *rdta;
  1428.   { /* getbasefreqs */
  1429.     int  ch;
  1430.  
  1431.     if (rdta->freqread) return TRUE;
  1432.  
  1433.     ch = getc(INFILE);
  1434.     if (! ((ch == 'F') || (ch == 'f')))  (void) ungetc(ch, INFILE);
  1435.  
  1436.     if (fscanf(INFILE, "%lf%lf%lf%lf",
  1437.                & rdta->freqa, & rdta->freqc,
  1438.                & rdta->freqg, & rdta->freqt) != 4 ||
  1439.         findch('\n') == EOF) {
  1440.       printf("ERROR: Problem reading user base frequencies\n");
  1441.       return  FALSE;
  1442.       }
  1443.  
  1444.     return TRUE;
  1445.   } /* getbasefreqs */
  1446.  
  1447.  
  1448. boolean getyspace (rdta)
  1449.     rawDNA   *rdta;
  1450.   { /* getyspace */
  1451.     long   size;
  1452.     int    i;
  1453.     yType *y0;
  1454.  
  1455.     if (! (rdta->y = (yType **) Malloc((rdta->numsp + 1) * sizeof(yType *)))) {
  1456.       printf("ERROR: Unable to obtain space for data array pointers\n");
  1457.       return  FALSE;
  1458.       }
  1459.  
  1460.     size = 4 * (rdta->sites / 4 + 1);
  1461.     if (! (y0 = (yType *) Malloc((rdta->numsp + 1) * size * sizeof(yType)))) {
  1462.       printf("ERROR: Unable to obtain space for data array\n");
  1463.       return  FALSE;
  1464.       }
  1465.  
  1466.     for (i = 0; i <= rdta->numsp; i++) {
  1467.       rdta->y[i] = y0;
  1468.       y0 += size;
  1469.       }
  1470.  
  1471.     return  TRUE;
  1472.   } /* getyspace */
  1473.  
  1474.  
  1475. boolean setupTree (tr, nsites)
  1476.     tree  *tr;
  1477.     int    nsites;
  1478.   { /* setupTree */
  1479.     nodeptr  p0, p, q;
  1480.     int      i, j, tips, inter;
  1481.  
  1482.     tips  = tr->mxtips;
  1483.     inter = tr->mxtips - 1;
  1484.  
  1485.     if (!(p0 = (nodeptr) Malloc((tips + 3*inter) * sizeof(node)))) {
  1486.       printf("ERROR: Unable to obtain sufficient tree memory\n");
  1487.       return  FALSE;
  1488.       }
  1489.  
  1490.     if (!(tr->nodep = (nodeptr *) Malloc((2*tr->mxtips) * sizeof(nodeptr)))) {
  1491.       printf("ERROR: Unable to obtain sufficient tree memory, too\n");
  1492.       return  FALSE;
  1493.       }
  1494.  
  1495.     tr->nodep[0] = (node *) NULL;    /* Use as 1-based array */
  1496.  
  1497.     for (i = 1; i <= tips; i++) {    /* Set-up tips */
  1498.       p = p0++;
  1499.       p->x      = (xarray *) NULL;
  1500.       p->tip    = (yType *) NULL;
  1501.       p->number = i;
  1502.       p->next   = p;
  1503.       p->back   = (node *) NULL;
  1504.  
  1505.       tr->nodep[i] = p;
  1506.       }
  1507.  
  1508.     for (i = tips + 1; i <= tips + inter; i++) { /* Internal nodes */
  1509.       q = (node *) NULL;
  1510.       for (j = 1; j <= 3; j++) {
  1511.         p = p0++;
  1512.         p->x      = (xarray *) NULL;
  1513.         p->tip    = (yType *) NULL;
  1514.         p->number = i;
  1515.         p->next   = q;
  1516.         p->back   = (node *) NULL;
  1517.         q = p;
  1518.         }
  1519.       p->next->next->next = p;
  1520.       tr->nodep[i] = p;
  1521.       }
  1522.  
  1523.     tr->likelihood  = unlikely;
  1524.     tr->start       = (node *) NULL;
  1525.     tr->outgrnode   = tr->nodep[tr->outgr];
  1526.     tr->ntips       = 0;
  1527.     tr->nextnode    = 0;
  1528.     tr->opt_level   = 0;
  1529.     tr->prelabeled  = TRUE;
  1530.     tr->smoothed    = FALSE;
  1531.     tr->log_f_valid = 0;
  1532.  
  1533.     tr->log_f = (double *) Malloc(nsites * sizeof(double));
  1534.     if (! tr->log_f) {
  1535.       printf("ERROR: Unable to obtain sufficient tree memory, trey\n");
  1536.       return  FALSE;
  1537.       }
  1538.  
  1539.     return TRUE;
  1540.   } /* setupTree */
  1541.  
  1542.  
  1543. void freeTreeNode (p)       /* Free tree node (sector) associated data */
  1544.     nodeptr  p;
  1545.   { /* freeTreeNode */
  1546.     if (p) {
  1547.       if (p->x) {
  1548.         if (p->x->a) Free(p->x->a);
  1549.         Free(p->x);
  1550.         }
  1551.       }
  1552.   } /* freeTreeNode */
  1553.  
  1554.  
  1555. void freeTree (tr)
  1556.     tree  *tr;
  1557.   { /* freeTree */
  1558.     nodeptr  p, q;
  1559.     int  i, tips, inter;
  1560.  
  1561.     tips  = tr->mxtips;
  1562.     inter = tr->mxtips - 1;
  1563.  
  1564.     for (i = 1; i <= tips; i++) freeTreeNode(tr->nodep[i]);
  1565.  
  1566.     for (i = tips + 1; i <= tips + inter; i++) {
  1567.       if (p = tr->nodep[i]) {
  1568.         if (q = p->next) {
  1569.           freeTreeNode(q->next);
  1570.           freeTreeNode(q);
  1571.           }
  1572.         freeTreeNode(p);
  1573.         }
  1574.       }
  1575.  
  1576.     Free(tr->nodep[1]);       /* Free the actual nodes */
  1577.   } /* freeTree */
  1578.  
  1579.  
  1580. boolean getdata (adef, rdta, tr)
  1581.     analdef  *adef;
  1582.     rawDNA   *rdta;
  1583.     tree     *tr;
  1584.     /* read sequences */
  1585.   { /* getdata */
  1586.     int   i, j, k, l, basesread, basesnew, ch;
  1587.     int   meaning[256];          /*  meaning of input characters */ 
  1588.     char *nameptr;
  1589.     boolean  allread, firstpass;
  1590.  
  1591.     for (i = 0; i <= 255; i++) meaning[i] = 0;
  1592.     meaning['A'] =  1;
  1593.     meaning['B'] = 14;
  1594.     meaning['C'] =  2;
  1595.     meaning['D'] = 13;
  1596.     meaning['G'] =  4;
  1597.     meaning['H'] = 11;
  1598.     meaning['K'] = 12;
  1599.     meaning['M'] =  3;
  1600.     meaning['N'] = 15;
  1601.     meaning['O'] = 15;
  1602.     meaning['R'] =  5;
  1603.     meaning['S'] =  6;
  1604.     meaning['T'] =  8;
  1605.     meaning['U'] =  8;
  1606.     meaning['V'] =  7;
  1607.     meaning['W'] =  9;
  1608.     meaning['X'] = 15;
  1609.     meaning['Y'] = 10;
  1610.     meaning['?'] = 15;
  1611.     meaning['-'] = 15;
  1612.  
  1613.     basesread = basesnew = 0;
  1614.  
  1615.     allread = FALSE;
  1616.     firstpass = TRUE;
  1617.     ch = ' ';
  1618.  
  1619.     while (! allread) {
  1620.       for (i = 1; i <= tr->mxtips; i++) {     /*  Read data line */
  1621.  
  1622.         if (firstpass) {                      /*  Read species names */
  1623.           j = 1;
  1624.           while (whitechar(ch = getc(INFILE))) {  /*  Skip blank lines */
  1625.             if (ch == '\n')  j = 1;  else  j++;
  1626.             }
  1627.  
  1628.           if (j > nmlngth) {
  1629.             printf("ERROR: Blank name for species %d; ", i);
  1630.             printf("check number of species,\n");
  1631.             printf("       number of sites, and interleave option.\n");
  1632.             return  FALSE;
  1633.             }
  1634.  
  1635.           nameptr = tr->nodep[i]->name;
  1636.           for (k = 1; k < j; k++)  *nameptr++ = ' ';
  1637.  
  1638.           while (ch != '\n' && ch != EOF) {
  1639.             if (whitechar(ch))  ch = ' ';
  1640.             *nameptr++ = ch;
  1641.             if (++j > nmlngth) break;
  1642.             ch = getc(INFILE);
  1643.             }
  1644.  
  1645.           while (*(--nameptr) == ' ') ;          /*  remove trailing blanks */
  1646.           *(++nameptr) = '\0';                   /*  add null termination */
  1647.  
  1648.           if (ch == EOF) {
  1649.             printf("ERROR: End-of-file in name of species %d\n", i);
  1650.             return  FALSE;
  1651.             }
  1652.           }    /* if (firstpass) */
  1653.  
  1654.         j = basesread;
  1655.         while ((j < rdta->sites)
  1656.                && ((ch = getc(INFILE)) != EOF)
  1657.                && ((! adef->interleaved) || (ch != '\n'))) {
  1658.           uppercase(& ch);
  1659.           if (meaning[ch] || ch == '.') {
  1660.             j++;
  1661.             if (ch == '.') {
  1662.               if (i != 1) ch = rdta->y[1][j];
  1663.               else {
  1664.                 printf("ERROR: Dot (.) found at site %d of sequence 1\n", j);
  1665.                 return  FALSE;
  1666.                 }
  1667.               }
  1668.             rdta->y[i][j] = ch;
  1669.             }
  1670.           else if (whitechar(ch) || digitchar(ch)) ;
  1671.           else {
  1672.             printf("ERROR: Bad base (%c) at site %d of sequence %d\n",
  1673.                     ch, j, i);
  1674.             return  FALSE;
  1675.             }
  1676.           }
  1677.  
  1678.         if (ch == EOF) {
  1679.           printf("ERROR: End-of-file at site %d of sequence %d\n", j, i);
  1680.           return  FALSE;
  1681.           }
  1682.  
  1683.         if (! firstpass && (j == basesread)) i--;        /* no data on line */
  1684.         else if (i == 1) basesnew = j;
  1685.         else if (j != basesnew) {
  1686.           printf("ERROR: Sequences out of alignment\n");
  1687.           printf("%d (instead of %d) residues read in sequence %d\n",
  1688.                   j - basesread, basesnew - basesread, i);
  1689.           return  FALSE;
  1690.           }
  1691.  
  1692.         while (ch != '\n' && ch != EOF) ch = getc(INFILE);  /* flush line */
  1693.         }                                                  /* next sequence */
  1694.       firstpass = FALSE;
  1695.       basesread = basesnew;
  1696.       allread = (basesread >= rdta->sites);
  1697.       }
  1698.  
  1699.     /*  Print listing of sequence alignment */
  1700.  
  1701.     if (adef->prdata) {
  1702.       j = nmlngth - 5 + ((rdta->sites + ((rdta->sites - 1)/10))/2);
  1703.       if (j < nmlngth - 1) j = nmlngth - 1;
  1704.       if (j > 37) j = 37;
  1705.       printf("Name"); for (i=1;i<=j;i++) putchar(' '); printf("Sequences\n");
  1706.       printf("----"); for (i=1;i<=j;i++) putchar(' '); printf("---------\n");
  1707.       putchar('\n');
  1708.  
  1709.       for (i = 1; i <= rdta->sites; i += 60) {
  1710.         l = i + 59;
  1711.         if (l > rdta->sites) l = rdta->sites;
  1712.  
  1713.         if (adef->userwgt) {
  1714.           printf("Weights   ");
  1715.           for (j = 11; j <= nmlngth+3; j++) putchar(' ');
  1716.           for (k = i; k <= l; k++) {
  1717.             putchar(itobase36(rdta->wgt[k]));
  1718.             if (((k % 10) == 0) && (k < l)) putchar(' ');
  1719.             }
  1720.           putchar('\n');
  1721.           }
  1722.  
  1723.         if (rdta->categs > 1) {
  1724.           printf("Categories");
  1725.           for (j = 11; j <= nmlngth+3; j++) putchar(' ');
  1726.           for (k = i; k <= l; k++) {
  1727.             putchar(itobase36(rdta->sitecat[k]));
  1728.             if (((k % 10) == 0) && (k < l)) putchar(' ');
  1729.             }
  1730.           putchar('\n');
  1731.           }
  1732.  
  1733.         for (j = 1; j <= tr->mxtips; j++) {
  1734.           nameptr = tr->nodep[j]->name;
  1735.           k = nmlngth+3;
  1736.           while (ch = *nameptr++) {putchar(ch); k--;}
  1737.           while (--k >= 0) putchar(' ');
  1738.  
  1739.           for (k = i; k <= l; k++) {
  1740.             ch = rdta->y[j][k];
  1741.             if ((j > 1) && (ch == rdta->y[1][k])) ch = '.';
  1742.             putchar(ch);
  1743.             if (((k % 10) == 0) && (k < l)) putchar(' ');
  1744.             }
  1745.           putchar('\n');
  1746.           }
  1747.         putchar('\n');
  1748.         }
  1749.       }
  1750.  
  1751.     for (j = 1; j <= tr->mxtips; j++)    /* Convert characters to meanings */
  1752.       for (i = 1; i <= rdta->sites; i++) {
  1753.         rdta->y[j][i] = meaning[rdta->y[j][i]];
  1754.         }
  1755.  
  1756.     return  TRUE;
  1757.   } /* getdata */
  1758.  
  1759.  
  1760. boolean  getntrees (adef)
  1761.     analdef  *adef;
  1762.   { /* getntrees */
  1763.  
  1764.     if (fscanf(INFILE, "%d", &(adef->numutrees)) != 1 || findch('\n') == EOF) {
  1765.       printf("ERROR: Problem reading number of user trees\n");
  1766.       return  FALSE;
  1767.       }
  1768.  
  1769.     if (adef->nkeep == 0) adef->nkeep = adef->numutrees;
  1770.  
  1771.     return  TRUE;
  1772.   } /* getntrees */
  1773.  
  1774.  
  1775. boolean getinput (adef, rdta, cdta, tr)
  1776.     analdef      *adef;
  1777.     rawDNA       *rdta;
  1778.     crunchedDNA  *cdta;
  1779.     tree         *tr;
  1780.   { /* getinput */
  1781.     if (! getnums(rdta))                       return FALSE;
  1782.     if (! getoptions(adef, rdta, cdta, tr))    return FALSE;
  1783.     if (! adef->empf && ! getbasefreqs(rdta))  return FALSE;
  1784.     if (! getyspace(rdta))                     return FALSE;
  1785.     if (! setupTree(tr, rdta->sites))          return FALSE;
  1786.     if (! getdata(adef, rdta, tr))             return FALSE;
  1787.     if (adef->usertree && ! getntrees(adef))   return FALSE;
  1788.  
  1789.     return TRUE;
  1790.   } /* getinput */
  1791.  
  1792.  
  1793. void makeboot (adef, rdta, cdta)
  1794.     analdef      *adef;
  1795.     rawDNA       *rdta;
  1796.     crunchedDNA  *cdta;
  1797.   { /* makeboot */
  1798.     int  i, j, nonzero;
  1799.     double  randum();
  1800.  
  1801.     nonzero = 0;
  1802.     for (i = 1; i <= rdta->sites; i++)  if (rdta->wgt[i] > 0) nonzero++;
  1803.  
  1804.     for (j = 1; j <= nonzero; j++) cdta->aliaswgt[j] = 0;
  1805.     for (j = 1; j <= nonzero; j++)
  1806.       cdta->aliaswgt[(int) (nonzero*randum(& adef->boot)) + 1]++;
  1807.  
  1808.     j = 0;
  1809.     cdta->wgtsum = 0;
  1810.     for (i = 1; i <= rdta->sites; i++) {
  1811.       if (rdta->wgt[i] > 0)
  1812.         cdta->wgtsum += (rdta->wgt2[i] = rdta->wgt[i] * cdta->aliaswgt[++j]);
  1813.       else
  1814.         rdta->wgt2[i] = 0;
  1815.       }
  1816.   } /* makeboot */
  1817.  
  1818.  
  1819. void sitesort (rdta, cdta)
  1820.     rawDNA       *rdta;
  1821.     crunchedDNA  *cdta;
  1822.     /* Shell sort keeping sites, weights in same order */
  1823.   { /* sitesort */
  1824.     int  gap, i, j, jj, jg, k;
  1825.     boolean  flip, tied;
  1826.  
  1827.     for (gap = rdta->sites / 2; gap > 0; gap /= 2) {
  1828.       for (i = gap + 1; i <= rdta->sites; i++) {
  1829.         j = i - gap;
  1830.  
  1831.         do {
  1832.           jj = cdta->alias[j];
  1833.           jg = cdta->alias[j+gap];
  1834.           flip = (rdta->sitecat[jj] >  rdta->sitecat[jg]);
  1835.           tied = (rdta->sitecat[jj] == rdta->sitecat[jg]);
  1836.           for (k = 1; (k <= rdta->numsp) && tied; k++) {
  1837.             flip = (rdta->y[k][jj] >  rdta->y[k][jg]);
  1838.             tied = (rdta->y[k][jj] == rdta->y[k][jg]);
  1839.             }
  1840.           if (flip) {
  1841.             cdta->alias[j]     = jg;
  1842.             cdta->alias[j+gap] = jj;
  1843.             j -= gap;
  1844.             }
  1845.           } while (flip && (j > 0));
  1846.  
  1847.         }  /* for (i ...   */
  1848.       }    /* for (gap ... */
  1849.   } /* sitesort */
  1850.  
  1851.  
  1852. void sitecombcrunch (rdta, cdta)
  1853.     rawDNA       *rdta;
  1854.     crunchedDNA  *cdta;
  1855.     /* combine sites that have identical patterns (and nonzero weight) */
  1856.   { /* sitecombcrunch */
  1857.     int  i, sitei, j, sitej, k;
  1858.     boolean  tied;
  1859.  
  1860.     i = 0;
  1861.     cdta->alias[0] = cdta->alias[1];
  1862.     cdta->aliaswgt[0] = 0;
  1863.  
  1864.     for (j = 1; j <= rdta->sites; j++) {
  1865.       sitei = cdta->alias[i];
  1866.       sitej = cdta->alias[j];
  1867.       tied = (rdta->sitecat[sitei] == rdta->sitecat[sitej]);
  1868.  
  1869.       for (k = 1; tied && (k <= rdta->numsp); k++)
  1870.         tied = (rdta->y[k][sitei] == rdta->y[k][sitej]);
  1871.  
  1872.       if (tied) {
  1873.         cdta->aliaswgt[i] += rdta->wgt2[sitej];
  1874.         }
  1875.       else {
  1876.         if (cdta->aliaswgt[i] > 0) i++;
  1877.         cdta->aliaswgt[i] = rdta->wgt2[sitej];
  1878.         cdta->alias[i] = sitej;
  1879.         }
  1880.       }
  1881.  
  1882.     cdta->endsite = i;
  1883.     if (cdta->aliaswgt[i] > 0) cdta->endsite++;
  1884.   } /* sitecombcrunch */
  1885.  
  1886.  
  1887. boolean makeweights (adef, rdta, cdta)
  1888.     analdef      *adef;
  1889.     rawDNA       *rdta;
  1890.     crunchedDNA  *cdta;
  1891.     /* make up weights vector to avoid duplicate computations */
  1892.   { /* makeweights */
  1893.     int  i;
  1894.  
  1895.     if (adef->boot)  makeboot(adef, rdta, cdta);
  1896.     else  for (i = 1; i <= rdta->sites; i++)  rdta->wgt2[i] = rdta->wgt[i];
  1897.  
  1898.     for (i = 1; i <= rdta->sites; i++)  cdta->alias[i] = i;
  1899.     sitesort(rdta, cdta);
  1900.     sitecombcrunch(rdta, cdta);
  1901.  
  1902.     printf("Total weight of positions in analysis = %d\n", cdta->wgtsum);
  1903.     printf("There are %d distinct data patterns (columns)\n\n", cdta->endsite);
  1904.  
  1905.     return TRUE;
  1906.   } /* makeweights */
  1907.  
  1908.  
  1909. boolean makevalues (rdta, cdta)
  1910.     rawDNA       *rdta;
  1911.     crunchedDNA  *cdta;
  1912.     /* set up fractional likelihoods at tips */
  1913.   { /* makevalues */
  1914.     double  temp, wtemp;
  1915.     int  i, j;
  1916.  
  1917.     for (i = 1; i <= rdta->numsp; i++) {    /* Pack and move tip data */
  1918.       for (j = 0; j < cdta->endsite; j++) {
  1919.         rdta->y[i-1][j] = rdta->y[i][cdta->alias[j]];
  1920.         }
  1921.       }
  1922.  
  1923.     for (j = 0; j < cdta->endsite; j++) {
  1924.       cdta->patcat[j] = i = rdta->sitecat[cdta->alias[j]];
  1925.       cdta->patrat[j] = temp = rdta->catrat[i];
  1926.       cdta->wr[j]  = wtemp = temp * cdta->aliaswgt[j];
  1927.       cdta->wr2[j] = temp * wtemp;
  1928.       }
  1929.  
  1930.     return TRUE;
  1931.   } /* makevalues */
  1932.  
  1933.  
  1934. boolean empiricalfreqs (rdta, cdta)
  1935.     rawDNA       *rdta;
  1936.     crunchedDNA  *cdta;
  1937.     /* Get empirical base frequencies from the data */
  1938.   { /* empiricalfreqs */
  1939.     double  sum, suma, sumc, sumg, sumt, wj, fa, fc, fg, ft;
  1940.     int     i, j, k, code;
  1941.     yType  *yptr;
  1942.  
  1943.     rdta->freqa = 0.25;
  1944.     rdta->freqc = 0.25;
  1945.     rdta->freqg = 0.25;
  1946.     rdta->freqt = 0.25;
  1947.     for (k = 1; k <= 8; k++) {
  1948.       suma = 0.0;
  1949.       sumc = 0.0;
  1950.       sumg = 0.0;
  1951.       sumt = 0.0;
  1952.       for (i = 0; i < rdta->numsp; i++) {
  1953.         yptr = rdta->y[i];
  1954.         for (j = 0; j < cdta->endsite; j++) {
  1955.           code = *yptr++;
  1956.           fa = rdta->freqa * ( code       & 1);
  1957.           fc = rdta->freqc * ((code >> 1) & 1);
  1958.           fg = rdta->freqg * ((code >> 2) & 1);
  1959.           ft = rdta->freqt * ((code >> 3) & 1);
  1960.           wj = cdta->aliaswgt[j] / (fa + fc + fg + ft);
  1961.           suma += wj * fa;
  1962.           sumc += wj * fc;
  1963.           sumg += wj * fg;
  1964.           sumt += wj * ft;
  1965.           }
  1966.         }
  1967.       sum = suma + sumc + sumg + sumt;
  1968.       rdta->freqa = suma / sum;
  1969.       rdta->freqc = sumc / sum;
  1970.       rdta->freqg = sumg / sum;
  1971.       rdta->freqt = sumt / sum;
  1972.       }
  1973.  
  1974.     return TRUE;
  1975.   } /* empiricalfreqs */
  1976.  
  1977.  
  1978. void reportfreqs (adef, rdta)
  1979.     analdef      *adef;
  1980.     rawDNA       *rdta;
  1981.   { /* reportfreqs */
  1982.     double  suma, sumb;
  1983.  
  1984.     if (adef->empf) printf("Empirical ");
  1985.     printf("Base Frequencies:\n\n");
  1986.     printf("   A    %10.5f\n",   rdta->freqa);
  1987.     printf("   C    %10.5f\n",   rdta->freqc);
  1988.     printf("   G    %10.5f\n",   rdta->freqg);
  1989.     printf("  T(U)  %10.5f\n\n", rdta->freqt);
  1990.     rdta->freqr = rdta->freqa + rdta->freqg;
  1991.     rdta->invfreqr = 1.0/rdta->freqr;
  1992.     rdta->freqar = rdta->freqa * rdta->invfreqr;
  1993.     rdta->freqgr = rdta->freqg * rdta->invfreqr;
  1994.     rdta->freqy = rdta->freqc + rdta->freqt;
  1995.     rdta->invfreqy = 1.0/rdta->freqy;
  1996.     rdta->freqcy = rdta->freqc * rdta->invfreqy;
  1997.     rdta->freqty = rdta->freqt * rdta->invfreqy;
  1998.     printf("Transition/transversion ratio = %10.6f\n\n", rdta->ttratio);
  1999.     suma = rdta->ttratio*rdta->freqr*rdta->freqy
  2000.          - (rdta->freqa*rdta->freqg + rdta->freqc*rdta->freqt);
  2001.     sumb = rdta->freqa*rdta->freqgr + rdta->freqc*rdta->freqty;
  2002.     rdta->xi = suma/(suma+sumb);
  2003.     rdta->xv = 1.0 - rdta->xi;
  2004.     if (rdta->xi <= 0.0) {
  2005.       printf("WARNING: This transition/transversion ratio\n");
  2006.       printf("         is impossible with these base frequencies!\n");
  2007.       printf("Transition/transversion parameter reset\n\n");
  2008.       rdta->xi = 0.000001;
  2009.       rdta->xv = 1.0 - rdta->xi;
  2010.       rdta->ttratio = (sumb * rdta->xi / rdta->xv 
  2011.                        + rdta->freqa * rdta->freqg
  2012.                        + rdta->freqc * rdta->freqt)
  2013.                     / (rdta->freqr * rdta->freqy);
  2014.       printf("Transition/transversion ratio = %10.6f\n\n", rdta->ttratio);
  2015.       }
  2016.     printf("(Transition/transversion parameter = %10.6f)\n\n",
  2017.             rdta->xi/rdta->xv);
  2018.     rdta->fracchange = 2.0 * rdta->xi * (rdta->freqa * rdta->freqgr
  2019.                                        + rdta->freqc * rdta->freqty)
  2020.                      + rdta->xv * (1.0 - rdta->freqa * rdta->freqa
  2021.                                        - rdta->freqc * rdta->freqc
  2022.                                        - rdta->freqg * rdta->freqg
  2023.                                        - rdta->freqt * rdta->freqt);
  2024.   } /* reportfreqs */
  2025.  
  2026.  
  2027. boolean linkdata2tree (rdta, cdta, tr)
  2028.     rawDNA       *rdta;
  2029.     crunchedDNA  *cdta;
  2030.     tree         *tr;
  2031.     /* Link data array to the tree tips */
  2032.   { /* linkdata2tree */
  2033.     int  i;
  2034.  
  2035.     for (i = 1; i <= tr->mxtips; i++) {    /* Associate data with tips */
  2036.       tr->nodep[i]->tip = &(rdta->y[i-1][0]);
  2037.       }
  2038.  
  2039.     tr->rdta       = rdta;
  2040.     tr->cdta       = cdta;
  2041.     return TRUE;
  2042.   } /* linkdata2tree */
  2043.  
  2044.  
  2045. xarray *setupxarray (npat)
  2046.     int npat;
  2047.   { /* setupxarray */
  2048.     xarray  *x;
  2049.     xtype   *data;
  2050.  
  2051.     x = (xarray *) Malloc(sizeof(xarray));
  2052.     if (x) {
  2053.       data = (xtype *) Malloc(4 * npat * sizeof(xtype));
  2054.       if (data) {
  2055.         x->a = data;
  2056.         x->c = data += npat;
  2057.         x->g = data += npat;
  2058.         x->t = data +  npat;
  2059.         x->prev = x->next = x;
  2060.         x->owner = (node *) NULL;
  2061.         }
  2062.       else {
  2063.         Free(x);
  2064.         return (xarray *) NULL;
  2065.         }
  2066.       }
  2067.     return x;
  2068.   } /* setupxarray */
  2069.  
  2070.  
  2071. boolean linkxarray (req, min, npat, freexptr, usedxptr)
  2072.     int  req, min, npat;
  2073.     xarray **freexptr, **usedxptr;
  2074.     /*  Link a set of xarrays */
  2075.   { /* linkxarray */
  2076.     xarray  *first, *prev, *x;
  2077.     int  i;
  2078.  
  2079.     first = prev = (xarray *) NULL;
  2080.     i = 0;
  2081.  
  2082.     do {
  2083.       x = setupxarray(npat);
  2084.       if (x) {
  2085.         if (! first) first = x;
  2086.         else {
  2087.           prev->next = x;
  2088.           x->prev = prev;
  2089.           }
  2090.         prev = x;
  2091.         i++;
  2092.         }
  2093.       else {
  2094.         printf("ERROR: Failure to get requested xarray memory\n");
  2095.         if (i < min)  return  FALSE;
  2096.         }
  2097.       } while ((i < req) && x);
  2098.  
  2099.     if (first) {
  2100.       first->prev = prev;
  2101.       prev->next = first;
  2102.       }
  2103.  
  2104.     *freexptr = first;
  2105.     *usedxptr = (xarray *) NULL;
  2106.  
  2107.     return  TRUE;
  2108.   } /* linkxarray */
  2109.  
  2110.  
  2111. boolean setupnodex (tr)
  2112.     tree  *tr;
  2113.   { /* setupnodex */
  2114.     nodeptr  p;
  2115.     int  i;
  2116.  
  2117.     for (i = tr->mxtips + 1; (i <= 2*(tr->mxtips) - 2); i++) {
  2118.       p = tr->nodep[i];
  2119.       if (! (p->x = setupxarray(tr->cdta->endsite))) {
  2120.         printf("ERROR: Failure to get internal node xarray memory\n");
  2121.         return  FALSE;
  2122.         }
  2123.       }
  2124.  
  2125.     return  TRUE;
  2126.   } /* setupnodex */
  2127.  
  2128.  
  2129. xarray *getxtip (p)
  2130.     nodeptr  p;
  2131.   { /* getxtip */
  2132.     xarray  *new;
  2133.     boolean  splice;
  2134.  
  2135.     if (! p) return (xarray *) NULL;
  2136.  
  2137.     splice = FALSE;
  2138.  
  2139.     if (p->x) {                  /* array is there; move to tail of list */
  2140.       new = p->x;
  2141.       if (new == new->prev) ;             /* linked to self; leave it */
  2142.       else if (new == usedxtip) usedxtip = usedxtip->next; /* at head */
  2143.       else if (new == usedxtip->prev) ;   /* already at tail */
  2144.       else {                              /* move to tail of list */
  2145.         new->prev->next = new->next;
  2146.         new->next->prev = new->prev;
  2147.         splice = TRUE;
  2148.         }
  2149.       }
  2150.  
  2151.     else if (freextip) {                 /* take from unused list */
  2152.       p->x = new = freextip;
  2153.       new->owner = p;
  2154.       if (new->prev != new) {            /* not only member of freelist */
  2155.         new->prev->next = new->next;
  2156.         new->next->prev = new->prev;
  2157.         freextip = new->next;
  2158.         }
  2159.       else
  2160.         freextip = (xarray *) NULL;
  2161.  
  2162.       splice = TRUE;
  2163.       }
  2164.  
  2165.     else if (usedxtip) {                 /* take from head of used list */
  2166.       usedxtip->owner->x = (xarray *) NULL;
  2167.       p->x = new = usedxtip;
  2168.       new->owner = p;
  2169.       usedxtip = usedxtip->next;
  2170.       }
  2171.  
  2172.     else {
  2173.       printf("ERROR: Unable to locate memory for tip %d.\n", p->number);
  2174.       return  (xarray *) NULL;
  2175.       }
  2176.  
  2177.     if (splice) {
  2178.       if (usedxtip) {                  /* list is not empty */
  2179.         usedxtip->prev->next = new;
  2180.         new->prev = usedxtip->prev;
  2181.         usedxtip->prev = new;
  2182.         new->next = usedxtip;
  2183.         }
  2184.       else
  2185.         usedxtip = new->prev = new->next = new;
  2186.       }
  2187.  
  2188.     return  new;
  2189.   } /* getxtip */
  2190.  
  2191.  
  2192. xarray *getxnode (p)
  2193.     nodeptr  p;
  2194.     /* Ensure that internal node p has memory */
  2195.   { /* getxnode */
  2196.     nodeptr  s;
  2197.  
  2198.     if (! (p->x)) {  /*  Move likelihood array on this node to sector p */
  2199.       if ((s = p->next)->x || (s = s->next)->x) {
  2200.         p->x = s->x;
  2201.         s->x = (xarray *) NULL;
  2202.         }
  2203.       else {
  2204.         printf("ERROR: Unable to locate memory at node %d.\n", p->number);
  2205.         exit(1);
  2206.         }
  2207.       }
  2208.     return  p->x;
  2209.   } /* getxnode */
  2210.  
  2211.  
  2212. boolean newview (tr, p)                      /*  Update likelihoods at node */
  2213.     tree    *tr;
  2214.     nodeptr  p;
  2215.   { /* newview */
  2216.     double   z1, lz1, xvlz1, z2, lz2, xvlz2;
  2217.     nodeptr  q, r;
  2218.     xtype   *x1a, *x1c, *x1g, *x1t, *x2a, *x2c, *x2g, *x2t,
  2219.             *x3a, *x3c, *x3g, *x3t;
  2220.     int  i;
  2221.  
  2222.     if (p->tip) {             /*  Make sure that data are at tip */
  2223.       int code;
  2224.       yType *yptr;
  2225.  
  2226.       if (p->x) return TRUE;  /*  They are already there */
  2227.  
  2228.       if (! getxtip(p)) return FALSE; /*  They are not, so get memory */
  2229.       x3a = &(p->x->a[0]);    /*  Move tip data to xarray */
  2230.       x3c = &(p->x->c[0]);
  2231.       x3g = &(p->x->g[0]);
  2232.       x3t = &(p->x->t[0]);
  2233.       yptr = p->tip;
  2234.       for (i = 0; i < tr->cdta->endsite; i++) {
  2235.         code = *yptr++;
  2236.         *x3a++ =  code       & 1;
  2237.         *x3c++ = (code >> 1) & 1;
  2238.         *x3g++ = (code >> 2) & 1;
  2239.         *x3t++ = (code >> 3) & 1;
  2240.         }
  2241.       return TRUE;
  2242.       }
  2243.  
  2244. /*  Internal node needs update */
  2245.  
  2246.     q = p->next->back;
  2247.     r = p->next->next->back;
  2248.  
  2249.     while ((! p->x) || (! q->x) || (! r->x)) {
  2250.       if (! q->x) if (! newview(tr, q))  return FALSE;
  2251.       if (! r->x) if (! newview(tr, r))  return FALSE;
  2252.       if (! p->x) if (! getxnode(p)) return FALSE;
  2253.       }
  2254.  
  2255.     x1a = &(q->x->a[0]);
  2256.     x1c = &(q->x->c[0]);
  2257.     x1g = &(q->x->g[0]);
  2258.     x1t = &(q->x->t[0]);
  2259.     z1 = q->z;
  2260.     lz1 = (z1 > zmin) ? log(z1) : log(zmin);
  2261.     xvlz1 = tr->rdta->xv * lz1;
  2262.  
  2263.     x2a = &(r->x->a[0]);
  2264.     x2c = &(r->x->c[0]);
  2265.     x2g = &(r->x->g[0]);
  2266.     x2t = &(r->x->t[0]);
  2267.     z2 = r->z;
  2268.     lz2 = (z2 > zmin) ? log(z2) : log(zmin);
  2269.     xvlz2 = tr->rdta->xv * lz2;
  2270.  
  2271.     x3a = &(p->x->a[0]);
  2272.     x3c = &(p->x->c[0]);
  2273.     x3g = &(p->x->g[0]);
  2274.     x3t = &(p->x->t[0]);
  2275.  
  2276.     { double  zz1table[maxcategories+1], zv1table[maxcategories+1],
  2277.               zz2table[maxcategories+1], zv2table[maxcategories+1],
  2278.               *zz1ptr, *zv1ptr, *zz2ptr, *zv2ptr, *rptr;
  2279.       double  fx1r, fx1y, fx1n, suma1, sumg1, sumc1, sumt1,
  2280.               fx2r, fx2y, fx2n, ki, tempi, tempj;
  2281.       int  *cptr;
  2282.  
  2283. #     ifdef Vectorize
  2284.         double zz1[maxpatterns], zv1[maxpatterns],
  2285.                zz2[maxpatterns], zv2[maxpatterns];
  2286.         int  cat;
  2287. #     else
  2288.         double zz1, zv1, zz2, zv2;
  2289.         int cat;
  2290. #     endif
  2291.  
  2292.       rptr   = &(tr->rdta->catrat[1]);
  2293.       zz1ptr = &(zz1table[1]);
  2294.       zv1ptr = &(zv1table[1]);
  2295.       zz2ptr = &(zz2table[1]);
  2296.       zv2ptr = &(zv2table[1]);
  2297.  
  2298. #     ifdef Vectorize
  2299. #       pragma IVDEP
  2300. #     endif
  2301.  
  2302.       for (i = 1; i <= tr->rdta->categs; i++) {   /* exps for each category */
  2303.         ki = *rptr++;
  2304.         *zz1ptr++ = exp(ki *   lz1);
  2305.         *zv1ptr++ = exp(ki * xvlz1);
  2306.         *zz2ptr++ = exp(ki *   lz2);
  2307.         *zv2ptr++ = exp(ki * xvlz2);
  2308.         }
  2309.  
  2310.       cptr = &(tr->cdta->patcat[0]);
  2311.  
  2312. #     ifdef Vectorize
  2313. #       pragma IVDEP
  2314.         for (i = 0; i < tr->cdta->endsite; i++) {
  2315.           cat = *cptr++;
  2316.           zz1[i] = zz1table[cat];
  2317.           zv1[i] = zv1table[cat];
  2318.           zz2[i] = zz2table[cat];
  2319.           zv2[i] = zv2table[cat];
  2320.         }
  2321.  
  2322. #       pragma IVDEP
  2323.         for (i = 0; i < tr->cdta->endsite; i++) {
  2324.           fx1r = tr->rdta->freqa * *x1a + tr->rdta->freqg * *x1g;
  2325.           fx1y = tr->rdta->freqc * *x1c + tr->rdta->freqt * *x1t;
  2326.           fx1n = fx1r + fx1y;
  2327.           tempi = fx1r * tr->rdta->invfreqr;
  2328.           tempj = zv1[i] * (tempi-fx1n) + fx1n;
  2329.           suma1 = zz1[i] * (*x1a++ - tempi) + tempj;
  2330.           sumg1 = zz1[i] * (*x1g++ - tempi) + tempj;
  2331.           tempi = fx1y * tr->rdta->invfreqy;
  2332.           tempj = zv1[i] * (tempi-fx1n) + fx1n;
  2333.           sumc1 = zz1[i] * (*x1c++ - tempi) + tempj;
  2334.           sumt1 = zz1[i] * (*x1t++ - tempi) + tempj;
  2335.  
  2336.           fx2r = tr->rdta->freqa * *x2a + tr->rdta->freqg * *x2g;
  2337.           fx2y = tr->rdta->freqc * *x2c + tr->rdta->freqt * *x2t;
  2338.           fx2n = fx2r + fx2y;
  2339.           tempi = fx2r * tr->rdta->invfreqr;
  2340.           tempj = zv2[i] * (tempi-fx2n) + fx2n;
  2341.           *x3a++ = suma1 * (zz2[i] * (*x2a++ - tempi) + tempj);
  2342.           *x3g++ = sumg1 * (zz2[i] * (*x2g++ - tempi) + tempj);
  2343.           tempi = fx2y * tr->rdta->invfreqy;
  2344.           tempj = zv2[i] * (tempi-fx2n) + fx2n;
  2345.           *x3c++ = sumc1 * (zz2[i] * (*x2c++ - tempi) + tempj);
  2346.           *x3t++ = sumt1 * (zz2[i] * (*x2t++ - tempi) + tempj);
  2347.           }
  2348.  
  2349. #     else  /*  Not Vectorize  */
  2350.         for (i = 0; i < tr->cdta->endsite; i++) {
  2351.           cat = *cptr++;
  2352.           zz1 = zz1table[cat];
  2353.           zv1 = zv1table[cat];
  2354.           fx1r = tr->rdta->freqa * *x1a + tr->rdta->freqg * *x1g;
  2355.           fx1y = tr->rdta->freqc * *x1c + tr->rdta->freqt * *x1t;
  2356.           fx1n = fx1r + fx1y;
  2357.           tempi = fx1r * tr->rdta->invfreqr;
  2358.           tempj = zv1 * (tempi-fx1n) + fx1n;
  2359.           suma1 = zz1 * (*x1a++ - tempi) + tempj;
  2360.           sumg1 = zz1 * (*x1g++ - tempi) + tempj;
  2361.           tempi = fx1y * tr->rdta->invfreqy;
  2362.           tempj = zv1 * (tempi-fx1n) + fx1n;
  2363.           sumc1 = zz1 * (*x1c++ - tempi) + tempj;
  2364.           sumt1 = zz1 * (*x1t++ - tempi) + tempj;
  2365.  
  2366.           zz2 = zz2table[cat];
  2367.           zv2 = zv2table[cat];
  2368.           fx2r = tr->rdta->freqa * *x2a + tr->rdta->freqg * *x2g;
  2369.           fx2y = tr->rdta->freqc * *x2c + tr->rdta->freqt * *x2t;
  2370.           fx2n = fx2r + fx2y;
  2371.           tempi = fx2r * tr->rdta->invfreqr;
  2372.           tempj = zv2 * (tempi-fx2n) + fx2n;
  2373.           *x3a++ = suma1 * (zz2 * (*x2a++ - tempi) + tempj);
  2374.           *x3g++ = sumg1 * (zz2 * (*x2g++ - tempi) + tempj);
  2375.           tempi = fx2y * tr->rdta->invfreqy;
  2376.           tempj = zv2 * (tempi-fx2n) + fx2n;
  2377.           *x3c++ = sumc1 * (zz2 * (*x2c++ - tempi) + tempj);
  2378.           *x3t++ = sumt1 * (zz2 * (*x2t++ - tempi) + tempj);
  2379.           }
  2380. #     endif  /*  Vectorize or not  */
  2381.  
  2382.       return TRUE;
  2383.       }
  2384.   } /* newview */
  2385.  
  2386.  
  2387. double evaluate (tr, p)
  2388.     tree    *tr;
  2389.     nodeptr  p;
  2390.   { /* evaluate */
  2391.     double   sum, z, lz, xvlz,
  2392.              ki, fx1a, fx1c, fx1g, fx1t, fx1r, fx1y, fx2r, fx2y,
  2393.              suma, sumb, sumc, term;
  2394.  
  2395. #   ifdef Vectorize
  2396.        double   zz[maxpatterns], zv[maxpatterns];
  2397. #   else
  2398.        double   zz,zv;
  2399. #   endif
  2400.  
  2401.     double   zztable[maxcategories+1], zvtable[maxcategories+1],
  2402.             *zzptr, *zvptr;
  2403.     double  *log_f, *rptr;
  2404.     xtype   *x1a, *x1c, *x1g, *x1t, *x2a, *x2c, *x2g, *x2t;
  2405.     nodeptr  q;
  2406.     int  cat, *cptr, i, *wptr;
  2407.  
  2408.     q = p->back;
  2409.     while ((! p->x) || (! q->x)) {
  2410.       if (! (p->x)) if (! newview(tr, p)) return badEval;
  2411.       if (! (q->x)) if (! newview(tr, q)) return badEval;
  2412.       }
  2413.  
  2414.     x1a = &(p->x->a[0]);
  2415.     x1c = &(p->x->c[0]);
  2416.     x1g = &(p->x->g[0]);
  2417.     x1t = &(p->x->t[0]);
  2418.  
  2419.     x2a = &(q->x->a[0]);
  2420.     x2c = &(q->x->c[0]);
  2421.     x2g = &(q->x->g[0]);
  2422.     x2t = &(q->x->t[0]);
  2423.  
  2424.     z = p->z;
  2425.     if (z < zmin) z = zmin;
  2426.     lz = log(z);
  2427.     xvlz = tr->rdta->xv * lz;
  2428.  
  2429.     rptr  = &(tr->rdta->catrat[1]);
  2430.     zzptr = &(zztable[1]);
  2431.     zvptr = &(zvtable[1]);
  2432.  
  2433. #   ifdef Vectorize
  2434. #     pragma IVDEP
  2435. #   endif
  2436.  
  2437.     for (i = 1; i <= tr->rdta->categs; i++) {
  2438.       ki = *rptr++;
  2439.       *zzptr++ = exp(ki *   lz);
  2440.       *zvptr++ = exp(ki * xvlz);
  2441.       }
  2442.  
  2443.     wptr = &(tr->cdta->aliaswgt[0]);
  2444.     cptr = &(tr->cdta->patcat[0]);
  2445.     log_f = tr->log_f;
  2446.     tr->log_f_valid = tr->cdta->endsite;
  2447.     sum = 0.0;
  2448.  
  2449. #   ifdef Vectorize
  2450. #     pragma IVDEP
  2451.       for (i = 0; i < tr->cdta->endsite; i++) {
  2452.         cat   = *cptr++;
  2453.         zz[i] = zztable[cat];
  2454.         zv[i] = zvtable[cat];
  2455.       }
  2456.  
  2457. #     pragma IVDEP
  2458.       for (i = 0; i < tr->cdta->endsite; i++) {
  2459.         fx1a = tr->rdta->freqa * *x1a++;
  2460.         fx1g = tr->rdta->freqg * *x1g++;
  2461.         fx1c = tr->rdta->freqc * *x1c++;
  2462.         fx1t = tr->rdta->freqt * *x1t++;
  2463.         suma = fx1a * *x2a + fx1c * *x2c + fx1g * *x2g + fx1t * *x2t;
  2464.         fx2r = tr->rdta->freqa * *x2a++ + tr->rdta->freqg * *x2g++;
  2465.         fx2y = tr->rdta->freqc * *x2c++ + tr->rdta->freqt * *x2t++;
  2466.         fx1r = fx1a + fx1g;
  2467.         fx1y = fx1c + fx1t;
  2468.         sumc = (fx1r + fx1y) * (fx2r + fx2y);
  2469.         sumb = fx1r * fx2r * tr->rdta->invfreqr + fx1y * fx2y * tr->rdta->invfreqy;
  2470.         suma -= sumb;
  2471.         sumb -= sumc;
  2472.         term = log(zz[i] * suma + zv[i] * sumb + sumc);
  2473.         sum += *wptr++ * term;
  2474.         *log_f++ = term;
  2475.         }
  2476.  
  2477. #   else  /*  Not Vectorize  */
  2478.       for (i = 0; i < tr->cdta->endsite; i++) {
  2479.         cat  = *cptr++;
  2480.         zz   = zztable[cat];
  2481.         zv   = zvtable[cat];
  2482.         fx1a = tr->rdta->freqa * *x1a++;
  2483.         fx1g = tr->rdta->freqg * *x1g++;
  2484.         fx1c = tr->rdta->freqc * *x1c++;
  2485.         fx1t = tr->rdta->freqt * *x1t++;
  2486.         suma = fx1a * *x2a + fx1c * *x2c + fx1g * *x2g + fx1t * *x2t;
  2487.         fx2r = tr->rdta->freqa * *x2a++ + tr->rdta->freqg * *x2g++;
  2488.         fx2y = tr->rdta->freqc * *x2c++ + tr->rdta->freqt * *x2t++;
  2489.         fx1r = fx1a + fx1g;
  2490.         fx1y = fx1c + fx1t;
  2491.         sumc = (fx1r + fx1y) * (fx2r + fx2y);
  2492.         sumb = fx1r * fx2r * tr->rdta->invfreqr + fx1y * fx2y * tr->rdta->invfreqy;
  2493.         suma -= sumb;
  2494.         sumb -= sumc;
  2495.         term = log(zz * suma + zv * sumb + sumc);
  2496.         sum += *wptr++ * term;
  2497.         *log_f++ = term;
  2498.         }
  2499. #   endif  /*  Vectorize or not  */
  2500.  
  2501.     tr->likelihood = sum;
  2502.     return  sum;
  2503.   } /* evaluate */
  2504.  
  2505.  
  2506.  
  2507.  
  2508. double makenewz (tr, p, q, z0, maxiter)
  2509.     tree    *tr;
  2510.     nodeptr  p, q;
  2511.     double   z0;
  2512.     int  maxiter;
  2513.   { /* makenewz */
  2514.     double  *abi, *bci, *sumci, *abptr, *bcptr, *sumcptr;
  2515.     double   dlnLidlz, dlnLdlz, d2lnLdlz2, z, zprev, zstep, lz, xvlz,
  2516.              ki, suma, sumb, sumc, ab, bc, inv_Li, t1, t2,
  2517.              fx1a, fx1c, fx1g, fx1t, fx1r, fx1y, fx2r, fx2y;
  2518.     double   zztable[maxcategories+1], zvtable[maxcategories+1];
  2519.     double   *zzptr, *zvptr;
  2520.     double  *rptr, *wrptr, *wr2ptr;
  2521.     xtype   *x1a, *x1c, *x1g, *x1t, *x2a, *x2c, *x2g, *x2t;
  2522.     int      cat, *cptr, i, curvatOK;
  2523.  
  2524.  
  2525.     while ((! p->x) || (! q->x)) {
  2526.       if (! (p->x)) if (! newview(tr, p)) return badZ;
  2527.       if (! (q->x)) if (! newview(tr, q)) return badZ;
  2528.       }
  2529.  
  2530.     x1a = &(p->x->a[0]);
  2531.     x1c = &(p->x->c[0]);
  2532.     x1g = &(p->x->g[0]);
  2533.     x1t = &(p->x->t[0]);
  2534.     x2a = &(q->x->a[0]);
  2535.     x2c = &(q->x->c[0]);
  2536.     x2g = &(q->x->g[0]);
  2537.     x2t = &(q->x->t[0]);
  2538.  
  2539.     { unsigned scratch_size;
  2540.       scratch_size = sizeof(double) * tr->cdta->endsite;
  2541.       if ((abi   = (double *) Malloc(scratch_size)) &&
  2542.           (bci   = (double *) Malloc(scratch_size)) &&
  2543.           (sumci = (double *) Malloc(scratch_size))) ;
  2544.       else {
  2545.         printf("ERROR: makenewz unable to obtain space for arrays\n");
  2546.         return badZ;
  2547.         }
  2548.       }
  2549.     abptr   = abi;
  2550.     bcptr   = bci;
  2551.     sumcptr = sumci;
  2552.  
  2553. #   ifdef Vectorize
  2554. #     pragma IVDEP
  2555. #   endif
  2556.  
  2557.     for (i = 0; i < tr->cdta->endsite; i++) {
  2558.       fx1a = tr->rdta->freqa * *x1a++;
  2559.       fx1g = tr->rdta->freqg * *x1g++;
  2560.       fx1c = tr->rdta->freqc * *x1c++;
  2561.       fx1t = tr->rdta->freqt * *x1t++;
  2562.       suma = fx1a * *x2a + fx1c * *x2c + fx1g * *x2g + fx1t * *x2t;
  2563.       fx2r = tr->rdta->freqa * *x2a++ + tr->rdta->freqg * *x2g++;
  2564.       fx2y = tr->rdta->freqc * *x2c++ + tr->rdta->freqt * *x2t++;
  2565.       fx1r = fx1a + fx1g;
  2566.       fx1y = fx1c + fx1t;
  2567.       *sumcptr++ = sumc = (fx1r + fx1y) * (fx2r + fx2y);
  2568.       sumb       = fx1r * fx2r * tr->rdta->invfreqr + fx1y * fx2y * tr->rdta->invfreqy;
  2569.       *abptr++   = suma - sumb;
  2570.       *bcptr++   = sumb - sumc;
  2571.       }
  2572.  
  2573.     z = z0;
  2574.     do {
  2575.       zprev = z;
  2576.       zstep = (1.0 - zmax) * z + zmin;
  2577.       curvatOK = FALSE;
  2578.  
  2579.       do {
  2580.         if (z < zmin) z = zmin;
  2581.         else if (z > zmax) z = zmax;
  2582.  
  2583.         lz    = log(z);
  2584.         xvlz  = tr->rdta->xv * lz;
  2585.         rptr  = &(tr->rdta->catrat[1]);
  2586.         zzptr = &(zztable[1]);
  2587.         zvptr = &(zvtable[1]);
  2588.  
  2589. #       ifdef Vectorize
  2590. #         pragma IVDEP
  2591. #       endif
  2592.  
  2593.         for (i = 1; i <= tr->rdta->categs; i++) {
  2594.           ki = *rptr++;
  2595.           *zzptr++ = exp(ki *   lz);
  2596.           *zvptr++ = exp(ki * xvlz);
  2597.           }
  2598.  
  2599.         abptr   = abi;
  2600.         bcptr   = bci;
  2601.         sumcptr = sumci;
  2602.         cptr    = &(tr->cdta->patcat[0]);
  2603.         wrptr   = &(tr->cdta->wr[0]);
  2604.         wr2ptr  = &(tr->cdta->wr2[0]);
  2605.         dlnLdlz = 0.0;                 /*  = d(ln(likelihood))/d(lz) */
  2606.         d2lnLdlz2 = 0.0;               /*  = d2(ln(likelihood))/d(lz)2 */
  2607.  
  2608. #       ifdef Vectorize
  2609. #         pragma IVDEP
  2610. #       endif
  2611.  
  2612.         for (i = 0; i < tr->cdta->endsite; i++) {
  2613.           cat    = *cptr++;              /*  ratecategory(i) */
  2614.           ab     = *abptr++ * zztable[cat];
  2615.           bc     = *bcptr++ * zvtable[cat];
  2616.           sumc   = *sumcptr++;
  2617.           inv_Li = 1.0/(ab + bc + sumc);
  2618.           t1     = ab * inv_Li;
  2619.           t2     = tr->rdta->xv * bc * inv_Li;
  2620.           dlnLidlz   = t1 + t2;
  2621.           dlnLdlz   += *wrptr++  * dlnLidlz;
  2622.           d2lnLdlz2 += *wr2ptr++ * (t1 + tr->rdta->xv * t2 - dlnLidlz * dlnLidlz);
  2623.           }
  2624.  
  2625.         if ((d2lnLdlz2 >= 0.0) && (z < zmax))
  2626.           zprev = z = 0.37 * z + 0.63;  /*  Bad curvature, shorten branch */
  2627.         else
  2628.           curvatOK = TRUE;
  2629.  
  2630.         } while (! curvatOK);
  2631.  
  2632.       if (d2lnLdlz2 < 0.0) {
  2633.         z *= exp(-dlnLdlz / d2lnLdlz2);
  2634.         if (z < zmin) z = zmin;
  2635.         if (z > 0.25 * zprev + 0.75)    /*  Limit steps toward z = 1.0 */
  2636.           z = 0.25 * zprev + 0.75;
  2637.         }
  2638.  
  2639.       if (z > zmax) z = zmax;
  2640.  
  2641.       } while ((--maxiter > 0) && (ABS(z - zprev) > zstep));
  2642.  
  2643.     Free(abi);
  2644.     Free(bci);
  2645.     Free(sumci);
  2646.  
  2647.     return  z;
  2648.   } /* makenewz */
  2649.  
  2650.  
  2651. boolean update (tr, p)
  2652.     tree    *tr;
  2653.     nodeptr  p;
  2654.   { /* update */
  2655.     nodeptr  q;
  2656.     double   z0, z;
  2657.  
  2658.     q = p->back;
  2659.     z0 = q->z;
  2660.     if ((z = makenewz(tr, p, q, z0, newzpercycle)) == badZ) return FALSE;
  2661.     p->z = q->z = z;
  2662.     if (ABS(z - z0) > deltaz)  tr->smoothed = FALSE;
  2663.     return TRUE;
  2664.   } /* update */
  2665.  
  2666.  
  2667. boolean smooth (tr, p)
  2668.     tree    *tr;
  2669.     nodeptr  p;
  2670.   { /* smooth */
  2671.     nodeptr  q;
  2672.  
  2673.     if (! update(tr, p))               return FALSE; /*  Adjust branch */
  2674.     if (! p->tip) {                                  /*  Adjust descendants */
  2675.         q = p->next;
  2676.         while (q != p) {
  2677.           if (! smooth(tr, q->back))   return FALSE;
  2678.           q = q->next;
  2679.           }
  2680.  
  2681. #     if ReturnSmoothedView
  2682.         if (! newview(tr, p)) return FALSE;
  2683. #     endif
  2684.       }
  2685.  
  2686.     return TRUE;
  2687.   } /* smooth */
  2688.  
  2689.  
  2690. boolean smoothTree (tr, maxtimes)
  2691.     tree  *tr;
  2692.     int    maxtimes;
  2693.   { /* smoothTree */
  2694.     nodeptr  p, q;
  2695.  
  2696.     p = tr->start;
  2697.  
  2698.     while (--maxtimes >= 0) {
  2699.       tr->smoothed = TRUE;
  2700.       if (! smooth(tr, p->back))       return FALSE;
  2701.       if (! p->tip) {
  2702.         q = p->next;
  2703.         while (q != p) {
  2704.           if (! smooth(tr, q->back))   return FALSE;
  2705.           q = q->next;
  2706.           }
  2707.         }
  2708.       if (tr->smoothed)  break;
  2709.       }
  2710.  
  2711.     return TRUE;
  2712.   } /* smoothTree */
  2713.  
  2714.  
  2715. boolean localSmooth (tr, p, maxtimes)    /* Smooth branches around p */
  2716.     tree    *tr;
  2717.     nodeptr  p;
  2718.     int  maxtimes;
  2719.   { /* localSmooth */
  2720.     nodeptr  q;
  2721.  
  2722.     if (p->tip) return FALSE;            /* Should be an error */
  2723.  
  2724.     while (--maxtimes >= 0) {
  2725.       tr->smoothed = TRUE;
  2726.       q = p;
  2727.       do {
  2728.         if (! update(tr, q)) return FALSE;
  2729.         q = q->next;
  2730.         } while (q != p);
  2731.       if (tr->smoothed)  break;
  2732.       }
  2733.  
  2734.     tr->smoothed = FALSE;             /* Only smooth locally */
  2735.     return TRUE;
  2736.   } /* localSmooth */
  2737.  
  2738.  
  2739. void hookup (p, q, z)
  2740.     nodeptr  p, q;
  2741.     double   z;
  2742.   { /* hookup */
  2743.     p->back = q;
  2744.     q->back = p;
  2745.     p->z = q->z = z;
  2746.   } /* hookup */
  2747.  
  2748.  
  2749. boolean insert (tr, p, q, glob)   /* Insert node p into branch q <-> q->back */
  2750.     tree    *tr;
  2751.     nodeptr  p, q;
  2752.     boolean  glob;             /* Smooth tree globally? */
  2753.  
  2754. /*                q
  2755.                  /.
  2756.              add/ .
  2757.                /  .
  2758.              pn   .
  2759.     s ---- p      .remove
  2760.              pnn  .
  2761.                \  .
  2762.              add\ .
  2763.                  \.      pn  = p->next;
  2764.                   r      pnn = p->next->next;
  2765.  */
  2766.  
  2767.   { /* insert */
  2768.     nodeptr  r, s;
  2769.  
  2770.     r = q->back;
  2771.     s = p->back;
  2772.  
  2773. #   if BestInsertAverage && ! Master
  2774.     { double  zqr, zqs, zrs, lzqr, lzqs, lzrs, lzsum, lzq, lzr, lzs, lzmax;
  2775.  
  2776.       if ((zqr = makenewz(tr, q, r, q->z,     iterations)) == badZ) return FALSE;
  2777.       if ((zqs = makenewz(tr, q, s, defaultz, iterations)) == badZ) return FALSE;
  2778.       if ((zrs = makenewz(tr, r, s, defaultz, iterations)) == badZ) return FALSE;
  2779.  
  2780.       lzqr = (zqr > zmin) ? log(zqr) : log(zmin);  /* long branches */
  2781.       lzqs = (zqs > zmin) ? log(zqs) : log(zmin);
  2782.       lzrs = (zrs > zmin) ? log(zrs) : log(zmin);
  2783.       lzsum = 0.5 * (lzqr + lzqs + lzrs);
  2784.  
  2785.       lzq = lzsum - lzrs;
  2786.       lzr = lzsum - lzqs;
  2787.       lzs = lzsum - lzqr;
  2788.       lzmax = log(zmax);
  2789.  
  2790.       if      (lzq > lzmax) {lzq = lzmax; lzr = lzqr; lzs = lzqs;} /* short */
  2791.       else if (lzr > lzmax) {lzr = lzmax; lzq = lzqr; lzs = lzrs;}
  2792.       else if (lzs > lzmax) {lzs = lzmax; lzq = lzqs; lzr = lzrs;}
  2793.  
  2794.       hookup(p->next,       q, exp(lzq));
  2795.       hookup(p->next->next, r, exp(lzr));
  2796.       hookup(p,             s, exp(lzs));
  2797.       }
  2798.  
  2799. #   else
  2800.     { double  z;
  2801.       z = sqrt(q->z);
  2802.       hookup(p->next,       q, z);
  2803.       hookup(p->next->next, r, z);
  2804.       }
  2805.  
  2806. #   endif
  2807.  
  2808.     if (! newview(tr, p)) return FALSE;   /*  So that p is valid at update */
  2809.     tr->opt_level = 0;
  2810.  
  2811. #   if ! Master         /*  Smoothings are done by slave */
  2812.       if (glob) {                                    /* Smooth whole tree */
  2813.         if (! smoothTree(tr, smoothings)) return FALSE;
  2814.         }
  2815.       else {                                         /* Smooth locale of p */
  2816.         if (! localSmooth(tr, p, smoothings)) return FALSE;
  2817.         }
  2818.  
  2819. #   else
  2820.       tr->likelihood = unlikely;
  2821. #   endif
  2822.     return  TRUE;
  2823.   } /* insert */
  2824.  
  2825.  
  2826. nodeptr  removeNode (tr, p)
  2827.     tree    *tr;
  2828.     nodeptr  p;
  2829.  
  2830. /*                q
  2831.                  .|
  2832.           remove. |
  2833.                .  |
  2834.              pn   |
  2835.     s ---- p      |add
  2836.              pnn  |
  2837.                .  |
  2838.           remove. |
  2839.                  .|      pn  = p->next;
  2840.                   r      pnn = p->next->next;
  2841.  */
  2842.  
  2843.     /* remove p and return where it was */
  2844.   { /* removeNode */
  2845.     double   zqr;
  2846.     nodeptr  q, r;
  2847.  
  2848.     q = p->next->back;
  2849.     r = p->next->next->back;
  2850.     zqr = q->z * r->z;
  2851. #   if ! Master
  2852.       if ((zqr = makenewz(tr, q, r, zqr, iterations)) == badZ) return (node *) NULL;
  2853. #   endif
  2854.     hookup(q, r, zqr);
  2855.  
  2856.     p->next->next->back = p->next->back = (node *) NULL;
  2857.     return  q;
  2858.   } /* removeNode */
  2859.  
  2860.  
  2861. boolean initrav (tr, p)
  2862.     tree    *tr;
  2863.     nodeptr  p;
  2864.   { /* initrav */
  2865.     nodeptr  q;
  2866.  
  2867.     if (! p->tip) {
  2868.       q = p->next;
  2869.       do {
  2870.         if (! initrav(tr, q->back))  return FALSE;
  2871.         q = q->next;
  2872.         } while (q != p);
  2873.       if (! newview(tr, p))         return FALSE;
  2874.       }
  2875.  
  2876.     return TRUE;
  2877.   } /* initrav */
  2878.  
  2879.  
  2880. nodeptr buildNewTip (tr, p)
  2881.     tree  *tr;
  2882.     nodeptr  p;
  2883.   { /* buildNewTip */
  2884.     nodeptr  q;
  2885.  
  2886.     q = tr->nodep[(tr->nextnode)++];
  2887.     hookup(p, q, defaultz);
  2888.     return  q;
  2889.   } /* buildNewTip */
  2890.  
  2891.  
  2892. boolean buildSimpleTree (tr, ip, iq, ir)
  2893.     tree  *tr;
  2894.     int    ip, iq, ir;
  2895.   { /* buildSimpleTree */
  2896.     /* p, q and r are tips meeting at s */
  2897.     nodeptr  p, s;
  2898.     int  i;
  2899.  
  2900.     i = MIN(ip, iq);
  2901.     if (ir < i)  i = ir; 
  2902.     tr->start = tr->nodep[i];
  2903.     tr->ntips = 3;
  2904.     p = tr->nodep[ip];
  2905.     hookup(p, tr->nodep[iq], defaultz);
  2906.     s = buildNewTip(tr, tr->nodep[ir]);
  2907.  
  2908.     return insert(tr, s, p, FALSE);  /* Smoothing is local to s */
  2909.   } /* buildSimpleTree */
  2910.  
  2911.  
  2912. #ifndef HasStrLib
  2913.  
  2914. char * strchr (str, chr)
  2915.     char *str;
  2916.     int   chr;
  2917.  { /* strchr */
  2918.     int  c;
  2919.  
  2920.     while (c = *str)  {if (c == chr) return str; str++;}
  2921.     return  (char *) NULL;
  2922.  } /* strchr */
  2923.  
  2924.  
  2925. char * strstr (str1, str2)
  2926.     char *str1, *str2;
  2927.  { /* strstr */
  2928.     char *s1, *s2;
  2929.     int  c;
  2930.  
  2931.     while (*(s1 = str1)) {
  2932.       s2 = str2;
  2933.       do {
  2934.         if (! (c = *s2++))  return str1;
  2935.         } 
  2936.         while (*s1++ == c);
  2937.       str1++;
  2938.       }
  2939.     return  (char *) NULL;
  2940.  } /* strstr */
  2941.  
  2942. #endif
  2943.  
  2944.  
  2945. boolean readKeyValue (string, key, format, value)
  2946.     char *string, *key, *format;
  2947.     void *value;
  2948.   { /* readKeyValue */
  2949.  
  2950.     if (! (string = (char*)strstr(string, key)))  return FALSE;
  2951.     string += strlen(key);
  2952.     if (! (string = (char*)strchr(string, '=')))  return FALSE;
  2953.     string++;
  2954.     return  sscanf(string, format, value);  /* 1 if read, otherwise 0 */
  2955.   } /* readKeyValue */
  2956.  
  2957.  
  2958. #if Master || Slave
  2959.  
  2960. double str_readTreeLikelihood (treestr)
  2961.     char *treestr;
  2962.   { /* str_readTreeLikelihood */
  2963.     double lk1;
  2964.     char    *com, *com_end;
  2965.     boolean  readKeyValue();
  2966.  
  2967.     if ((com = strchr(treestr, '[')) && (com < strchr(treestr, '('))
  2968.                                      && (com_end = strchr(com, ']'))) {
  2969.       com++;
  2970.       *com_end = 0;
  2971.       if (readKeyValue(com, likelihood_key, "%lg", (void *) &(lk1))) {
  2972.         *com_end = ']';
  2973.         return lk1;
  2974.         }
  2975.       }
  2976.  
  2977.     fprintf(stderr, "ERROR reading likelihood in receiveTree\n");
  2978.     return  badEval;
  2979.   } /* str_readTreeLikelihood */
  2980.  
  2981.  
  2982. boolean sendTree (comm, tr)
  2983.     comm_block  *comm;
  2984.     tree        *tr;
  2985.   { /* sendTree */
  2986.     char  *treestr;
  2987.     char  *treeString();
  2988. #   if Master
  2989.       void sendTreeNum();
  2990. #   endif
  2991.  
  2992.     comm->done_flag = tr->likelihood > 0.0;
  2993.     if (comm->done_flag)
  2994.       write_comm_msg(comm, NULL);
  2995.  
  2996.     else {
  2997.       treestr = (char *) Malloc((tr->ntips * (nmlngth+32)) + 256);
  2998.       if (! treestr) {
  2999.         fprintf(stderr, "sendTree: Malloc failure\n");
  3000.         return 0;
  3001.         }
  3002.  
  3003. #     if Master
  3004.         if (send_ahead >= MAX_SEND_AHEAD) {
  3005.           double new_likelihood;
  3006.           int  n_to_get;
  3007.  
  3008.           n_to_get = (send_ahead+1)/2;
  3009.           sendTreeNum(n_to_get);
  3010.           send_ahead -= n_to_get;
  3011.           read_comm_msg(& comm_slave, treestr);
  3012.           new_likelihood = str_readTreeLikelihood(treestr);
  3013.           if (new_likelihood == badEval)  return FALSE;
  3014.           if (! best_tr_recv || (new_likelihood > best_lk_recv)) {
  3015.             if (best_tr_recv)  Free(best_tr_recv);
  3016.             best_tr_recv = Malloc(strlen(treestr) + 1);
  3017.             strcpy(best_tr_recv, treestr);
  3018.             best_lk_recv = new_likelihood;
  3019.             }
  3020.           }
  3021.         send_ahead++;
  3022. #     endif           /*  End #if Master  */
  3023.  
  3024.       REPORT_SEND_TREE;
  3025.       (void) treeString(treestr, tr, tr->start->back, 1);
  3026.       write_comm_msg(comm, treestr);
  3027.  
  3028.       Free(treestr);
  3029.       }
  3030.  
  3031.     return TRUE;
  3032.   } /* sendTree */
  3033.  
  3034.  
  3035. boolean  receiveTree (comm, tr)
  3036.     comm_block  *comm;
  3037.     tree        *tr;
  3038.   { /* receiveTree */
  3039.     char   *treestr;
  3040.     boolean status;
  3041.     boolean str_treeReadLen();
  3042.  
  3043.     treestr = (char *) Malloc((tr->ntips * (nmlngth+32)) + 256);
  3044.     if (! treestr) {
  3045.       fprintf(stderr, "receiveTree: Malloc failure\n");
  3046.       return 0;
  3047.       }
  3048.  
  3049.     read_comm_msg(comm, treestr);
  3050.     if (comm->done_flag) {
  3051.       tr->likelihood = 1.0;
  3052.       status = TRUE;
  3053.       }
  3054.  
  3055.     else {
  3056. #     if Master
  3057.         if (best_tr_recv) {
  3058.           if (str_readTreeLikelihood(treestr) < best_lk_recv) {
  3059.             strcpy(treestr, best_tr_recv);  /* Overwrite new tree with best */
  3060.             }
  3061.           Free(best_tr_recv);
  3062.           best_tr_recv = NULL;
  3063.           }
  3064. #     endif           /*  End #if Master  */
  3065.  
  3066.       status = str_treeReadLen(treestr, tr);
  3067.       }
  3068.  
  3069.     Free(treestr);
  3070.     return status;
  3071.   } /* receiveTree */
  3072.  
  3073.  
  3074. void requestForWork ()
  3075.   { /* requestForWork */
  3076.     p4_send(DNAML_REQUEST, DNAML_DISPATCHER_ID, NULL, 0);
  3077.   } /* requestForWork */
  3078. #endif                  /* End #if Master || Slave  */
  3079.  
  3080.  
  3081. #if Master
  3082. void sendTreeNum(n_to_get)
  3083.     int n_to_get;
  3084.   { /* sendTreeNum */
  3085.     char scr[512];
  3086.  
  3087.     sprintf(scr, "%d", n_to_get);
  3088.     p4_send(DNAML_NUM_TREE, DNAML_MERGER_ID, scr, strlen(scr)+1);
  3089.   } /* sendTreeNum */
  3090.  
  3091.  
  3092. boolean  getReturnedTrees (tr, bt, n_tree_sent)
  3093.     tree     *tr;
  3094.     bestlist *bt;
  3095.     int n_tree_sent; /* number of trees sent to slaves */
  3096.   { /* getReturnedTrees */
  3097.     void sendTreeNum();
  3098.     boolean receiveTree();
  3099.  
  3100.     sendTreeNum(send_ahead);
  3101.     send_ahead = 0;
  3102.  
  3103.     if (! receiveTree(& comm_slave, tr))  return FALSE;
  3104.     tr->smoothed = TRUE;
  3105.     (void) saveBestTree(bt, tr);
  3106.  
  3107.     return TRUE;
  3108.   } /* getReturnedTrees */
  3109. #endif
  3110.  
  3111.  
  3112. void cacheZ (tr)
  3113.     tree  *tr;
  3114.   { /* cacheZ */
  3115.     nodeptr  p;
  3116.     int  nodes;
  3117.  
  3118.     nodes = tr->mxtips  +  3 * (tr->mxtips - 2);
  3119.     p = tr->nodep[1];
  3120.     while (nodes-- > 0) {p->z0 = p->z; p++;}
  3121.   } /* cacheZ */
  3122.  
  3123.  
  3124. void restoreZ (tr)
  3125.     tree  *tr;
  3126.   { /* restoreZ */
  3127.     nodeptr  p;
  3128.     int  nodes;
  3129.  
  3130.     nodes = tr->mxtips  +  3 * (tr->mxtips - 2);
  3131.     p = tr->nodep[1];
  3132.     while (nodes-- > 0) {p->z = p->z0; p++;}
  3133.   } /* restoreZ */
  3134.  
  3135.  
  3136. boolean testInsert (tr, p, q, bt, fast)
  3137.     tree     *tr;
  3138.     nodeptr   p, q;
  3139.     bestlist *bt;
  3140.     boolean   fast;
  3141.   { /* testInsert */
  3142.     double  qz;
  3143.     nodeptr  r;
  3144.  
  3145.     r = q->back;             /* Save original connection */
  3146.     qz = q->z;
  3147.     if (! insert(tr, p, q, ! fast)) return FALSE;
  3148.  
  3149. #   if ! Master
  3150.       if (evaluate(tr, fast ? p->next->next : tr->start) == badEval) return FALSE;
  3151.       (void) saveBestTree(bt, tr);
  3152. #   else  /* Master */
  3153.       tr->likelihood = unlikely;
  3154.       if (! sendTree(& comm_slave, tr))  return FALSE;
  3155. #   endif
  3156.  
  3157.     /* remove p from this branch */
  3158.  
  3159.     hookup(q, r, qz);
  3160.     p->next->next->back = p->next->back = (nodeptr) NULL;
  3161.     if (! fast) {            /* With fast add, other values are still OK */
  3162.       restoreZ(tr);          /*   Restore branch lengths */
  3163. #     if ! Master            /*   Regenerate x values */
  3164.         if (! initrav(tr, p->back))  return FALSE;
  3165.         if (! initrav(tr, q))        return FALSE;
  3166.         if (! initrav(tr, r))        return FALSE;
  3167. #     endif
  3168.       }
  3169.  
  3170.     return TRUE;
  3171.   } /* testInsert */
  3172.  
  3173.  
  3174. int addTraverse (tr, p, q, mintrav, maxtrav, bt, fast)
  3175.     tree     *tr;
  3176.     nodeptr   p, q;
  3177.     int       mintrav, maxtrav;
  3178.     bestlist *bt;
  3179.     boolean   fast;
  3180.   { /* addTraverse */
  3181.     int  tested, newtested;
  3182.  
  3183.     tested = 0;
  3184.     if (--mintrav <= 0) {           /* Moved minimum distance? */
  3185.       if (! testInsert(tr, p, q, bt, fast))  return badRear;
  3186.       tested++;
  3187.       }
  3188.  
  3189.     if ((! q->tip) && (--maxtrav > 0)) {    /* Continue traverse? */
  3190.       newtested = addTraverse(tr, p, q->next->back,
  3191.                               mintrav, maxtrav, bt, fast);
  3192.       if (newtested == badRear) return badRear;
  3193.       tested += newtested;
  3194.       newtested = addTraverse(tr, p, q->next->next->back,
  3195.                               mintrav, maxtrav, bt, fast);
  3196.       if (newtested == badRear) return badRear;
  3197.       tested += newtested;
  3198.       }
  3199.  
  3200.     return tested;
  3201.   } /* addTraverse */
  3202.  
  3203.  
  3204. int  rearrange (tr, p, mintrav, maxtrav, bt)
  3205.     tree     *tr;
  3206.     nodeptr   p;
  3207.     int       mintrav, maxtrav;
  3208.     bestlist *bt;
  3209.     /* rearranges the tree, globally or locally */
  3210.   { /* rearrange */
  3211.     double   p1z, p2z, q1z, q2z;
  3212.     nodeptr  p1, p2, q, q1, q2;
  3213.     int      tested, mintrav2, newtested;
  3214.  
  3215.     tested = 0;
  3216.     if (maxtrav < 1 || mintrav > maxtrav)  return tested;
  3217.  
  3218. /* Moving subtree forward in tree. */
  3219.  
  3220.     if (! p->tip) {
  3221.       p1 = p->next->back;
  3222.       p2 = p->next->next->back;
  3223.       if (! p1->tip || ! p2->tip) {
  3224.         p1z = p1->z;
  3225.         p2z = p2->z;
  3226.         if (! removeNode(tr, p)) return badRear;
  3227.         cacheZ(tr);
  3228.         if (! p1->tip) {
  3229.           newtested = addTraverse(tr, p, p1->next->back,
  3230.                                   mintrav, maxtrav, bt, FALSE);
  3231.           if (newtested == badRear) return badRear;
  3232.           tested += newtested;
  3233.           newtested = addTraverse(tr, p, p1->next->next->back,
  3234.                                   mintrav, maxtrav, bt, FALSE);
  3235.           if (newtested == badRear) return badRear;
  3236.           tested += newtested;
  3237.           }
  3238.  
  3239.         if (! p2->tip) {
  3240.           newtested = addTraverse(tr, p, p2->next->back,
  3241.                                   mintrav, maxtrav, bt, FALSE);
  3242.           if (newtested == badRear) return badRear;
  3243.           tested += newtested;
  3244.           newtested = addTraverse(tr, p, p2->next->next->back,
  3245.                                   mintrav, maxtrav, bt, FALSE);
  3246.           if (newtested == badRear) return badRear;
  3247.           tested += newtested;
  3248.           }
  3249.  
  3250.         hookup(p->next,       p1, p1z);  /*  Restore original tree */
  3251.         hookup(p->next->next, p2, p2z);
  3252.         if (! (initrav(tr, tr->start)
  3253.             && initrav(tr, tr->start->back))) return badRear;
  3254.         }
  3255.       }   /* if (! p->tip) */
  3256.  
  3257. /* Moving subtree backward in tree.  Minimum move is 2 to avoid duplicates */
  3258.  
  3259.     q = p->back;
  3260.     if (! q->tip && maxtrav > 1) {
  3261.       q1 = q->next->back;
  3262.       q2 = q->next->next->back;
  3263.       if (! q1->tip && (!q1->next->back->tip || !q1->next->next->back->tip) ||
  3264.           ! q2->tip && (!q2->next->back->tip || !q2->next->next->back->tip)) {
  3265.         q1z = q1->z;
  3266.         q2z = q2->z;
  3267.         if (! removeNode(tr, q)) return badRear;
  3268.         cacheZ(tr);
  3269.         mintrav2 = mintrav > 2 ? mintrav : 2;
  3270.  
  3271.         if (! q1->tip) {
  3272.           newtested = addTraverse(tr, q, q1->next->back,
  3273.                                   mintrav2 , maxtrav, bt, FALSE);
  3274.           if (newtested == badRear) return badRear;
  3275.           tested += newtested;
  3276.           newtested = addTraverse(tr, q, q1->next->next->back,
  3277.                                   mintrav2 , maxtrav, bt, FALSE);
  3278.           if (newtested == badRear) return badRear;
  3279.           tested += newtested;
  3280.           }
  3281.  
  3282.         if (! q2->tip) {
  3283.           newtested = addTraverse(tr, q, q2->next->back,
  3284.                                   mintrav2 , maxtrav, bt, FALSE);
  3285.           if (newtested == badRear) return badRear;
  3286.           tested += newtested;
  3287.           newtested = addTraverse(tr, q, q2->next->next->back,
  3288.                                   mintrav2 , maxtrav, bt, FALSE);
  3289.           if (newtested == badRear) return badRear;
  3290.           tested += newtested;
  3291.           }
  3292.  
  3293.         hookup(q->next,       q1, q1z);  /*  Restore original tree */
  3294.         hookup(q->next->next, q2, q2z);
  3295.         if (! (initrav(tr, tr->start)
  3296.             && initrav(tr, tr->start->back))) return badRear;
  3297.         }
  3298.       }   /* if (! q->tip && maxtrav > 1) */
  3299.  
  3300. /* Move other subtrees */
  3301.  
  3302.     if (! p->tip) {
  3303.       newtested = rearrange(tr, p->next->back,       mintrav, maxtrav, bt);
  3304.       if (newtested == badRear) return badRear;
  3305.       tested += newtested;
  3306.       newtested = rearrange(tr, p->next->next->back, mintrav, maxtrav, bt);
  3307.       if (newtested == badRear) return badRear;
  3308.       tested += newtested;
  3309.       }
  3310.  
  3311.     return  tested;
  3312.   } /* rearrange */
  3313.  
  3314. #ifdef NoGetPID 
  3315. /* macintosh */
  3316. /* this is NOT define we want, but what Mac define fits all compilers !? */
  3317. long getpid(void) 
  3318. {
  3319.     return 0;
  3320. }
  3321. #endif
  3322.  
  3323. FILE *fopen_pid (filenm, mode, name_pid)
  3324.     char *filenm, *mode, *name_pid;
  3325.   { /* fopen_pid */
  3326.  
  3327.     (void) sprintf(name_pid, "%s.%d", filenm, getpid());
  3328.     return  fopen(name_pid, mode);
  3329.   } /* fopen_pid */
  3330.  
  3331.  
  3332. #if DeleteCheckpointFile
  3333. void  unlink_pid (filenm)
  3334.     char *filenm;
  3335.   { /* unlink_pid */
  3336.     char scr[512];
  3337.  
  3338.     (void) sprintf(scr, "%s.%d", filenm, getpid());
  3339.     unlink(scr);
  3340.   } /* unlink_pid */
  3341. #endif
  3342.  
  3343.  
  3344. void  writeCheckpoint (tr)
  3345.     tree  *tr;
  3346.   { /* writeCheckpoint */
  3347.     char   filename[128];
  3348.     FILE  *checkpointf;
  3349.     void   treeOut();
  3350.  
  3351.     checkpointf = fopen_pid(checkpointname, "a", filename);
  3352.     if (checkpointf) {
  3353.       treeOut(checkpointf, tr, treeNewick);
  3354.       (void) fclose(checkpointf);
  3355.       }
  3356.   } /* writeCheckpoint */
  3357.  
  3358.  
  3359. node * findAnyTip(p)
  3360.     nodeptr  p;
  3361.   { /* findAnyTip */
  3362.     return  p->tip ? p : findAnyTip(p->next->back);
  3363.   } /* findAnyTip */
  3364.  
  3365.  
  3366. boolean  optimize (tr, maxtrav, bt)
  3367.     tree     *tr;
  3368.     int       maxtrav;
  3369.     bestlist *bt;
  3370.   { /* optimize */
  3371.     nodeptr  p;
  3372.     int    mintrav, tested;
  3373.  
  3374.     if (tr->ntips < 4)  return  TRUE;
  3375.  
  3376.     writeCheckpoint(tr);                    /* checkpoint the starting tree */
  3377.  
  3378.     if (maxtrav > tr->ntips - 3)  maxtrav = tr->ntips - 3;
  3379.     if (maxtrav <= tr->opt_level)  return  TRUE;
  3380.  
  3381.     printf("      Doing %s rearrangements\n",
  3382.              (maxtrav == 1)            ? "local" :
  3383.              (maxtrav < tr->ntips - 3) ? "regional" : "global");
  3384.  
  3385.     /* loop while tree gets better  */
  3386.  
  3387.     do {
  3388.       (void) startOpt(bt, tr);
  3389.       mintrav = tr->opt_level + 1;
  3390.  
  3391.       /* rearrange must start from a tip or it will miss some trees */
  3392.  
  3393.       p = findAnyTip(tr->start);
  3394.       tested = rearrange(tr, p->back, mintrav, maxtrav, bt);
  3395.       if (tested == badRear)  return FALSE;
  3396.  
  3397. #     if Master
  3398.         if (! getReturnedTrees(tr, bt, tested)) return FALSE;
  3399. #     endif
  3400.  
  3401.       bt->numtrees += tested;
  3402.       (void) setOptLevel(bt, maxtrav);
  3403.       if (! recallBestTree(bt, 1, tr)) return FALSE;   /* recover best tree */
  3404.  
  3405.       printf("      Tested %d alternative trees\n", tested);
  3406.       if (bt->improved) {
  3407.         printf("      Ln Likelihood =%14.5f\n", tr->likelihood);
  3408.         }
  3409.  
  3410.       writeCheckpoint(tr);                  /* checkpoint the new tree */
  3411.       } while (maxtrav > tr->opt_level);
  3412.  
  3413.     return TRUE;
  3414.   } /* optimize */
  3415.  
  3416.  
  3417. void coordinates (tr, p, lengthsum, tdptr)
  3418.     tree     *tr;
  3419.     nodeptr   p;
  3420.     double    lengthsum;
  3421.     drawdata *tdptr;
  3422.   { /* coordinates */
  3423.     /* establishes coordinates of nodes */
  3424.     double  x, z;
  3425.     nodeptr  q, first, last;
  3426.  
  3427.     if (p->tip) {
  3428.       p->xcoord = NINT(over * lengthsum);
  3429.       p->ymax = p->ymin = p->ycoord = tdptr->tipy;
  3430.       tdptr->tipy += down;
  3431.       if (lengthsum > tdptr->tipmax) tdptr->tipmax = lengthsum;
  3432.       }
  3433.  
  3434.     else {
  3435.       q = p->next;
  3436.       do {
  3437.         z = q->z;
  3438.         if (z < zmin) z = zmin;
  3439.         x = lengthsum - tr->rdta->fracchange * log(z);
  3440.         coordinates(tr, q->back, x, tdptr);
  3441.         q = q->next;
  3442.         } while (p == tr->start->back ? q != p->next : q != p);
  3443.  
  3444.       first = p->next->back;
  3445.       q = p;
  3446.       while (q->next != p) q = q->next;
  3447.       last = q->back;
  3448.       p->xcoord = NINT(over * lengthsum);
  3449.       p->ycoord = (first->ycoord + last->ycoord)/2;
  3450.       p->ymin = first->ymin;
  3451.       p->ymax = last->ymax;
  3452.       }
  3453.   } /* coordinates */
  3454.  
  3455.  
  3456. void drawline (tr, i, scale)
  3457.     tree   *tr;
  3458.     int     i;
  3459.     double  scale;
  3460.     /* draws one row of the tree diagram by moving up tree */
  3461.     /* Modified to handle 1000 taxa, October 16, 1991 */
  3462.   { /* drawline */
  3463.     nodeptr  p, q, r, first, last;
  3464.     int  n, j, k, l, extra;
  3465.     boolean  done;
  3466.  
  3467.     p = q = tr->start->back;
  3468.     extra = 0;
  3469.  
  3470.     if (i == p->ycoord) {
  3471.       k = q->number - tr->mxtips;
  3472.       for (j = k; j < 1000; j *= 10) putchar('-');
  3473.       printf("%d", k);
  3474.       extra = 1;
  3475.       }
  3476.     else printf("   ");
  3477.  
  3478.     do {
  3479.       if (! p->tip) {
  3480.         r = p->next;
  3481.         done = FALSE;
  3482.         do {
  3483.           if ((i >= r->back->ymin) && (i <= r->back->ymax)) {
  3484.             q = r->back;
  3485.             done = TRUE;
  3486.             }
  3487.           r = r->next;
  3488.           } while (! done && (p == tr->start->back ? r != p->next : r != p));
  3489.  
  3490.         first = p->next->back;
  3491.         r = p;
  3492.         while (r->next != p) r = r->next;
  3493.         last = r->back;
  3494.         if (p == tr->start->back) last = p->back;
  3495.         }
  3496.  
  3497.       done = (p->tip) || (p == q);
  3498.       n = NINT(scale*(q->xcoord - p->xcoord));
  3499.       if ((n < 3) && (! q->tip)) n = 3;
  3500.       n -= extra;
  3501.       extra = 0;
  3502.  
  3503.       if ((q->ycoord == i) && (! done)) {
  3504.         if (p->ycoord != q->ycoord) putchar('+');
  3505.         else                        putchar('-');
  3506.  
  3507.         if (! q->tip) {
  3508.           k = q->number - tr->mxtips;
  3509.           l = n - 3;
  3510.           for (j = k; j < 100; j *= 10)  l++;
  3511.           for (j = 1; j <= l; j++) putchar('-');
  3512.           printf("%d", k);
  3513.           extra = 1;
  3514.           }
  3515.         else for (j = 1; j <= n-1; j++) putchar('-');
  3516.         }
  3517.  
  3518.       else if (! p->tip) {
  3519.         if ((last->ycoord > i) && (first->ycoord < i) && (i != p->ycoord)) {
  3520.           putchar('!');
  3521.           for (j = 1; j <= n-1; j++) putchar(' ');
  3522.           }
  3523.         else for (j = 1; j <= n; j++) putchar(' ');
  3524.         }
  3525.  
  3526.       else
  3527.         for (j = 1; j <= n; j++) putchar(' ');
  3528.  
  3529.       p = q;
  3530.       } while (! done);
  3531.  
  3532.     if ((p->ycoord == i) && p->tip) {
  3533.       printf(" %s", p->name);
  3534.       }
  3535.  
  3536.     putchar('\n');
  3537.   } /* drawline */
  3538.  
  3539.  
  3540. void printTree (tr, adef)
  3541.     tree     *tr;
  3542.     analdef  *adef;
  3543.     /* prints out diagram of the tree */
  3544.   { /* printTree */
  3545.     drawdata  tipdata;
  3546.     double  scale;
  3547.     int  i, imax;
  3548.  
  3549.     if (adef->trprint) {
  3550.       putchar('\n');
  3551.       tipdata.tipy = 1;
  3552.       tipdata.tipmax = 0.0;
  3553.       coordinates(tr, tr->start->back, (double) 0.0, & tipdata);
  3554.       scale = 1.0 / tipdata.tipmax;
  3555.       imax = tipdata.tipy - down;
  3556.       for (i = 1; i <= imax; i++)  drawline(tr, i, scale);
  3557.       printf("\nRemember: ");
  3558.       if (adef->root) printf("(although rooted by outgroup) ");
  3559.       printf("this is an unrooted tree!\n\n");
  3560.       }
  3561.   } /* printTree */
  3562.  
  3563.  
  3564. double sigma (tr, p, sumlrptr)
  3565.     tree    *tr;
  3566.     nodeptr  p;
  3567.     double  *sumlrptr;
  3568.     /* compute standard deviation */
  3569.   { /* sigma */
  3570.     double  slope, sum, sumlr, z, zv, zz, lz,
  3571.             rat, suma, sumb, sumc, d2, d, li, temp, abzz, bczv, t3,
  3572.             fx1a, fx1c, fx1g, fx1t, fx1r, fx1y, fx2r, fx2y, w;
  3573.     double  *rptr;
  3574.     xtype   *x1a, *x1c, *x1g, *x1t, *x2a, *x2c, *x2g, *x2t;
  3575.     nodeptr  q;
  3576.     int  i, *wptr;
  3577.  
  3578.     q = p->back;
  3579.     while ((! p->x) || (! q->x)) {
  3580.       if (! (p->x)) if (! newview(tr, p)) return -1.0;
  3581.       if (! (q->x)) if (! newview(tr, q)) return -1.0;
  3582.       }
  3583.  
  3584.     x1a = &(p->x->a[0]);
  3585.     x1c = &(p->x->c[0]);
  3586.     x1g = &(p->x->g[0]);
  3587.     x1t = &(p->x->t[0]);
  3588.  
  3589.     x2a = &(q->x->a[0]);
  3590.     x2c = &(q->x->c[0]);
  3591.     x2g = &(q->x->g[0]);
  3592.     x2t = &(q->x->t[0]);
  3593.  
  3594.     z = p->z;
  3595.     if (z < zmin) z = zmin;
  3596.     lz = log(z);
  3597.  
  3598.     wptr = &(tr->cdta->aliaswgt[0]);
  3599.     rptr = &(tr->cdta->patrat[0]);
  3600.     sum = sumlr = slope = 0.0;
  3601.  
  3602. #   ifdef Vectorize
  3603. #     pragma IVDEP
  3604. #   endif
  3605.  
  3606.     for (i = 0; i < tr->cdta->endsite; i++) {
  3607.       rat  = *rptr++;
  3608.       zz   = exp(rat            * lz);
  3609.       zv   = exp(rat * tr->rdta->xv * lz);
  3610.  
  3611.       fx1a = tr->rdta->freqa * *x1a++;
  3612.       fx1g = tr->rdta->freqg * *x1g++;
  3613.       fx1c = tr->rdta->freqc * *x1c++;
  3614.       fx1t = tr->rdta->freqt * *x1t++;
  3615.       fx1r = fx1a + fx1g;
  3616.       fx1y = fx1c + fx1t;
  3617.       suma = fx1a * *x2a + fx1c * *x2c + fx1g * *x2g + fx1t * *x2t;
  3618.       fx2r = tr->rdta->freqa * *x2a++ + tr->rdta->freqg * *x2g++;
  3619.       fx2y = tr->rdta->freqc * *x2c++ + tr->rdta->freqt * *x2t++;
  3620.       sumc = (fx1r + fx1y) * (fx2r + fx2y);
  3621.       sumb = fx1r * fx2r * tr->rdta->invfreqr + fx1y * fx2y * tr->rdta->invfreqy;
  3622.       abzz = zz * (suma - sumb);
  3623.       bczv = zv * (sumb - sumc);
  3624.       li = sumc + abzz + bczv;
  3625.       t3 = tr->rdta->xv * bczv;
  3626.       d  = abzz + t3;
  3627.       d2 = rat * (abzz*(rat-1.0) + t3*(rat * tr->rdta->xv - 1.0));
  3628.  
  3629.       temp = rat * d / li;
  3630.       w = *wptr++;
  3631.       slope += w *  temp;
  3632.       sum   += w * (temp * temp - d2/li);
  3633.       sumlr += w * log(li / (suma + 1.0E-300));
  3634.       }
  3635.  
  3636.     *sumlrptr = sumlr;
  3637.     return (sum > 1.0E-300) ? z*(-slope + sqrt(slope*slope + 3.841*sum))/sum
  3638.                             : 1.0;
  3639.   } /* sigma */
  3640.  
  3641.  
  3642. void describe (tr, p)
  3643.     tree    *tr;
  3644.     nodeptr  p;
  3645.     /* print out information for one branch */
  3646.   { /* describe */
  3647.     double   z, s, sumlr;
  3648.     nodeptr  q;
  3649.     char    *nameptr;
  3650.     int      k, ch;
  3651.  
  3652.     q = p->back;
  3653.     printf("%4d          ", q->number - tr->mxtips);
  3654.     if (p->tip) {
  3655.       nameptr = p->name;
  3656.       k = nmlngth;
  3657.       while (ch = *nameptr++) {putchar(ch); k--;}
  3658.       while (--k >= 0) putchar(' ');
  3659.       }
  3660.     else {
  3661.       printf("%4d", p->number - tr->mxtips);
  3662.       for (k = 4; k < nmlngth; k++) putchar(' ');
  3663.       }
  3664.  
  3665.     z = q->z;
  3666.     if (z <= zmin) printf("    infinity");
  3667.     else printf("%15.5f", -log(z) * tr->rdta->fracchange);
  3668.  
  3669.     s = sigma(tr, q, & sumlr);
  3670.     printf("     (");
  3671.     if (z + s >= zmax) printf("     zero");
  3672.     else printf("%9.5f", (double) -log(z + s) * tr->rdta->fracchange);
  3673.     putchar(',');
  3674.     if (z - s <= zmin) printf("    infinity");
  3675.     else printf("%12.5f", (double) -log(z - s) * tr->rdta->fracchange);
  3676.     putchar(')');
  3677.  
  3678.     if      (sumlr > 2.995 ) printf(" **");
  3679.     else if (sumlr > 1.9205) printf(" *");
  3680.     putchar('\n');
  3681.  
  3682.     if (! p->tip) {
  3683.       describe(tr, p->next->back);
  3684.       describe(tr, p->next->next->back);
  3685.       }
  3686.   } /* describe */
  3687.  
  3688.  
  3689. void summarize (tr)
  3690.     tree  *tr;
  3691.     /* print out branch length information and node numbers */
  3692.   { /* summarize */
  3693.     printf("Ln Likelihood =%14.5f\n", tr->likelihood);
  3694.     putchar('\n');
  3695.     printf(" Between        And             Length");
  3696.     printf("      Approx. Confidence Limits\n");
  3697.     printf(" -------        ---             ------");
  3698.     printf("      ------- ---------- ------\n");
  3699.  
  3700.     describe(tr, tr->start->back->next->back);
  3701.     describe(tr, tr->start->back->next->next->back);
  3702.     describe(tr, tr->start);
  3703.     putchar('\n');
  3704.     printf("     *  = significantly positive, P < 0.05\n");
  3705.     printf("     ** = significantly positive, P < 0.01\n\n\n");
  3706.   } /* summarize */
  3707.  
  3708.  
  3709. /*===========  This is a problem if tr->start->back is a tip!  ===========*/
  3710. /*  All routines should be contrived so that tr->start->back is not a tip */
  3711.  
  3712. char *treeString (treestr, tr, p, form)
  3713.     char  *treestr;
  3714.     tree  *tr;
  3715.     nodeptr  p;
  3716.     int  form;
  3717.     /* write string with representation of tree */
  3718.     /*   form == 1 -> Newick tree */
  3719.     /*   form == 2 -> Prolog fact */
  3720.     /*   form == 3 -> PHYLIP tree */
  3721.   { /* treeString */
  3722.     double  x, z;
  3723.     char  *nameptr;
  3724.     int    c;
  3725.  
  3726.     if (p == tr->start->back) {
  3727.       if (form != treePHYLIP) {
  3728.         if (form == treeProlog) {
  3729.           (void) sprintf(treestr, "phylip_tree(");
  3730.           while (*treestr) treestr++;            /* move pointer to null */
  3731.           }
  3732.  
  3733.         (void) sprintf(treestr, "[&&%s: version = '%s'",
  3734.                                  programName, programVersion);
  3735.         while (*treestr) treestr++;
  3736.  
  3737.         (void) sprintf(treestr, ", %s = %15.13g",
  3738.                                  likelihood_key, tr->likelihood);
  3739.         while (*treestr) treestr++;
  3740.  
  3741.         (void) sprintf(treestr, ", %s = %d", ntaxa_key, tr->ntips);
  3742.         while (*treestr) treestr++;
  3743.  
  3744.         (void) sprintf(treestr,", %s = %d", opt_level_key, tr->opt_level);
  3745.         while (*treestr) treestr++;
  3746.  
  3747.         (void) sprintf(treestr, ", %s = %d", smoothed_key, tr->smoothed);
  3748.         while (*treestr) treestr++;
  3749.  
  3750.         (void) sprintf(treestr, "]%s", form == treeProlog ? ", " : " ");
  3751.         while (*treestr) treestr++;
  3752.         }
  3753.       }
  3754.  
  3755.     if (p->tip) {
  3756.       if (form != treePHYLIP) *treestr++ = '\'';
  3757.       nameptr = p->name;
  3758.       while (c = *nameptr++) {
  3759.         if (form != treePHYLIP) {if (c == '\'') *treestr++ = '\'';}
  3760.         else if (c == ' ') {c = '_';}
  3761.         *treestr++ = c;
  3762.         }
  3763.       if (form != treePHYLIP) *treestr++ = '\'';
  3764.       }
  3765.  
  3766.     else {
  3767.       *treestr++ = '(';
  3768.       treestr = treeString(treestr, tr, p->next->back, form);
  3769.       *treestr++ = ',';
  3770.       treestr = treeString(treestr, tr, p->next->next->back, form);
  3771.       if (p == tr->start->back) {
  3772.         *treestr++ = ',';
  3773.         treestr = treeString(treestr, tr, p->back, form);
  3774.         }
  3775.       *treestr++ = ')';
  3776.       }
  3777.  
  3778.     if (p == tr->start->back) {
  3779.       (void) sprintf(treestr, ":0.0%s\n", (form != treeProlog) ? ";" : ").");
  3780.       }
  3781.     else {
  3782.       z = p->z;
  3783.       if (z < zmin) z = zmin;
  3784.       x = -log(z) * tr->rdta->fracchange;
  3785.       (void) sprintf(treestr, ": %8.6f", x);  /* prolog needs the space */
  3786.       }
  3787.  
  3788.     while (*treestr) treestr++;     /* move pointer up to null termination */
  3789.     return  treestr;
  3790.   } /* treeString */
  3791.  
  3792.  
  3793. void treeOut (treefile, tr, form)
  3794.     FILE  *treefile;
  3795.     tree  *tr;
  3796.     int    form;
  3797.     /* write out file with representation of final tree */
  3798.   { /* treeOut */
  3799.     int    c;
  3800.     char  *cptr, *treestr;
  3801.  
  3802.     treestr = (char *) Malloc((tr->ntips * (nmlngth+32)) + 256);
  3803.     if (! treestr) {
  3804.       fprintf(stderr, "treeOut: Malloc failure\n");
  3805.       exit(1);
  3806.       }
  3807.  
  3808.     (void) treeString(treestr, tr, tr->start->back, form);
  3809.     cptr = treestr;
  3810.     while (c = *cptr++) putc(c, treefile);
  3811.  
  3812.     Free(treestr);
  3813.   } /* treeOut */
  3814.  
  3815.  
  3816. /*=======================================================================*/
  3817. /*                         Read a tree from a file                       */
  3818. /*=======================================================================*/
  3819.  
  3820.  
  3821. /*  1.0.A  Processing of quotation marks in comment removed
  3822.  */
  3823.  
  3824. int treeFinishCom (fp, strp)
  3825.     FILE   *fp;
  3826.     char  **strp;
  3827.   { /* treeFinishCom */
  3828.     int  ch;
  3829.  
  3830.     while ((ch = getc(fp)) != EOF && ch != ']') {
  3831.       if (strp != NULL) *(*strp)++ = ch;    /* save character  */
  3832.       if (ch == '[') {                      /* nested comment; find its end */
  3833.         if ((ch = treeFinishCom(fp, strp)) == EOF)  break;
  3834.         if (strp != NULL) *(*strp)++ = ch;  /* save closing ]  */
  3835.         }
  3836.       }
  3837.  
  3838.     if (strp != NULL) **strp = '\0';        /* terminate string  */
  3839.     return  ch;
  3840.   } /* treeFinishCom */
  3841.  
  3842.  
  3843. int treeGetCh (fp)
  3844.     FILE  *fp;                   /* get next nonblank, noncomment character */
  3845.   { /* treeGetCh */
  3846.     int  ch;
  3847.  
  3848.     while ((ch = getc(fp)) != EOF) {
  3849.       if (whitechar(ch)) ;
  3850.       else if (ch == '[') {                   /* comment; find its end */
  3851.         if ((ch = treeFinishCom(fp, (char **) NULL)) == EOF)  break;
  3852.         }
  3853.       else  break;
  3854.       }
  3855.  
  3856.     return  ch;
  3857.   } /* treeGetCh */
  3858.  
  3859.  
  3860. boolean  treeLabelEnd (ch)
  3861.       int ch;
  3862.   { /* treeLabelEnd */
  3863.     switch (ch) {
  3864.         case EOF:  case '\0':  case '\t':  case '\n':  case ' ':
  3865.         case ':':  case ',':   case '(':   case ')':   case '[':
  3866.         case ';':
  3867.           return TRUE;
  3868.         default:
  3869.           break;
  3870.         }
  3871.     return FALSE;
  3872.   } /* treeLabelEnd */
  3873.  
  3874.  
  3875. boolean  treeGetLabel (fp, lblPtr, maxlen)
  3876.     FILE  *fp;
  3877.     char  *lblPtr;
  3878.     int    maxlen;
  3879.   { /* treeGetLabel */
  3880.     int      ch;
  3881.     boolean  done, quoted, lblfound;
  3882.  
  3883.     if (--maxlen < 0) lblPtr = (char *) NULL;  /* reserves space for '\0' */
  3884.     else if (lblPtr == NULL) maxlen = 0;
  3885.  
  3886.     ch = getc(fp);
  3887.     done = treeLabelEnd(ch);
  3888.  
  3889.     lblfound = ! done;
  3890.     quoted = (ch == '\'');
  3891.     if (quoted && ! done) {ch = getc(fp); done = (ch == EOF);}
  3892.  
  3893.     while (! done) {
  3894.       if (quoted) {
  3895.         if (ch == '\'') {ch = getc(fp); if (ch != '\'') break;}
  3896.         }
  3897.  
  3898.       else if (treeLabelEnd(ch)) break;
  3899.  
  3900.       else if (ch == '_') ch = ' ';  /* unquoted _ goes to space */
  3901.  
  3902.       if (--maxlen >= 0) *lblPtr++ = ch;
  3903.       ch = getc(fp);
  3904.       if (ch == EOF) break;
  3905.       }
  3906.  
  3907.     if (ch != EOF)  (void) ungetc(ch, fp);
  3908.  
  3909.     if (lblPtr != NULL) *lblPtr = '\0';
  3910.  
  3911.     return lblfound;
  3912.   } /* treeGetLabel */
  3913.  
  3914.  
  3915. boolean  treeFlushLabel (fp)
  3916.     FILE  *fp;
  3917.   { /* treeFlushLabel */
  3918.     return  treeGetLabel(fp, (char *) NULL, (int) 0);
  3919.   } /* treeFlushLabel */
  3920.  
  3921.  
  3922. int  treeFindTipByLabel (str, tr)
  3923.     char  *str;  /*  label string pointer */
  3924.     tree  *tr;
  3925.   { /* treeFindTipByLabel */
  3926.     nodeptr  q;
  3927.     char    *nameptr;
  3928.     int      ch, i, n;
  3929.     boolean  found;
  3930.  
  3931.     for (n = 1; n <= tr->mxtips; n++) {
  3932.       q = tr->nodep[n];
  3933.       if (! (q->back)) {          /*  Only consider unused tips */
  3934.         i = 0;
  3935.         nameptr = q->name;
  3936.         while ((found = (str[i++] == (ch = *nameptr++))) && ch) ;
  3937.         if (found) return n;
  3938.         }
  3939.       }
  3940.  
  3941.     printf("ERROR: Cannot find tree species: %s\n", str);
  3942.  
  3943.     return  0;
  3944.   } /* treeFindTipByLabel */
  3945.  
  3946.  
  3947. int  treeFindTipName (fp, tr)
  3948.     FILE  *fp;
  3949.     tree  *tr;
  3950.   { /* treeFindTipName */
  3951.     char    *nameptr, str[nmlngth+2];
  3952.     int      n;
  3953.  
  3954.     if (tr->prelabeled) {
  3955.       if (treeGetLabel(fp, str, nmlngth+2))
  3956.         n = treeFindTipByLabel(str, tr);
  3957.       else
  3958.         n = 0;
  3959.       }
  3960.  
  3961.     else if (tr->ntips < tr->mxtips) {
  3962.       n = tr->ntips + 1;
  3963.       nameptr = tr->nodep[n]->name;
  3964.       if (! treeGetLabel(fp, nameptr, nmlngth+1)) n = 0;
  3965.       }
  3966.  
  3967.     else {
  3968.       n = 0;
  3969.       }
  3970.  
  3971.     return  n;
  3972.   } /* treeFindTipName */
  3973.  
  3974.  
  3975. void  treeEchoContext (fp1, fp2, n)
  3976.     FILE  *fp1, *fp2;
  3977.     int    n;
  3978.  { /* treeEchoContext */
  3979.    int      ch;
  3980.    boolean  waswhite;
  3981.  
  3982.    waswhite = TRUE;
  3983.  
  3984.    while (n > 0 && ((ch = getc(fp1)) != EOF)) {
  3985.      if (whitechar(ch)) {
  3986.        ch = waswhite ? '\0' : ' ';
  3987.        waswhite = TRUE;
  3988.        }
  3989.      else {
  3990.        waswhite = FALSE;
  3991.        }
  3992.  
  3993.      if (ch > '\0') {putc(ch, fp2); n--;}
  3994.      }
  3995.  } /* treeEchoContext */
  3996.  
  3997.  
  3998. boolean treeProcessLength (fp, dptr)
  3999.     FILE    *fp;
  4000.     double  *dptr;
  4001.   { /* treeProcessLength */
  4002.     int  ch;
  4003.  
  4004.     if ((ch = treeGetCh(fp)) == EOF)  return FALSE;    /*  Skip comments */
  4005.     (void) ungetc(ch, fp);
  4006.  
  4007.     if (fscanf(fp, "%lf", dptr) != 1) {
  4008.       printf("ERROR: treeProcessLength: Problem reading branch length\n");
  4009.       treeEchoContext(fp, stdout, 40);
  4010.       printf("\n");
  4011.       return  FALSE;
  4012.       }
  4013.  
  4014.     return  TRUE;
  4015.   } /* treeProcessLength */
  4016.  
  4017.  
  4018. boolean  treeFlushLen (fp)
  4019.     FILE  *fp;
  4020.   { /* treeFlushLen */
  4021.     double  dummy;
  4022.     int     ch;
  4023.  
  4024.     if ((ch = treeGetCh(fp)) == ':') return treeProcessLength(fp, & dummy);
  4025.  
  4026.     if (ch != EOF) (void) ungetc(ch, fp);
  4027.     return TRUE;
  4028.   } /* treeFlushLen */
  4029.  
  4030.  
  4031. boolean  treeNeedCh (fp, c1, where)
  4032.     FILE  *fp;
  4033.     int    c1;
  4034.     char  *where;
  4035.   { /* treeNeedCh */
  4036.     int  c2;
  4037.  
  4038.     if ((c2 = treeGetCh(fp)) == c1)  return TRUE;
  4039.  
  4040.     printf("ERROR: Expecting '%c' %s tree; found:", c1, where);
  4041.     if (c2 == EOF) {
  4042.       printf("End-of-File");
  4043.       }
  4044.     else {
  4045.       ungetc(c2, fp);
  4046.       treeEchoContext(fp, stdout, 40);
  4047.       }
  4048.     putchar('\n');
  4049.     return FALSE;
  4050.   } /* treeNeedCh */
  4051.  
  4052.  
  4053. boolean  addElementLen (fp, tr, p)
  4054.     FILE    *fp;
  4055.     tree    *tr;
  4056.     nodeptr  p;
  4057.   { /* addElementLen */
  4058.     double   z, branch;
  4059.     nodeptr  q;
  4060.     int      n, ch;
  4061.  
  4062.     if ((ch = treeGetCh(fp)) == '(') {     /*  A new internal node */
  4063.       n = (tr->nextnode)++;
  4064.       if (n > 2*(tr->mxtips) - 2) {
  4065.         if (tr->rooted || n > 2*(tr->mxtips) - 1) {
  4066.           printf("ERROR: Too many internal nodes.  Is tree rooted?\n");
  4067.           printf("       Deepest splitting should be a trifurcation.\n");
  4068.           return FALSE;
  4069.           }
  4070.         else {
  4071.           tr->rooted = TRUE;
  4072.           }
  4073.         }
  4074.       q = tr->nodep[n];
  4075.       if (! addElementLen(fp, tr, q->next))        return FALSE;
  4076.       if (! treeNeedCh(fp, ',', "in"))             return FALSE;
  4077.       if (! addElementLen(fp, tr, q->next->next))  return FALSE;
  4078.       if (! treeNeedCh(fp, ')', "in"))             return FALSE;
  4079.       (void) treeFlushLabel(fp);
  4080.       }
  4081.  
  4082.     else {                               /*  A new tip */
  4083.       ungetc(ch, fp);
  4084.       if ((n = treeFindTipName(fp, tr)) <= 0)          return FALSE;
  4085.       q = tr->nodep[n];
  4086.       if (tr->start->number > n)  tr->start = q;
  4087.       (tr->ntips)++;
  4088.       }                                  /* End of tip processing */
  4089.  
  4090.     if (tr->userlen) {
  4091.       if (! treeNeedCh(fp, ':', "in"))             return FALSE;
  4092.       if (! treeProcessLength(fp, & branch))       return FALSE;
  4093.       z = exp(-branch / tr->rdta->fracchange);
  4094.       if (z > zmax)  z = zmax;
  4095.       hookup(p, q, z);
  4096.       }
  4097.     else {
  4098.       if (! treeFlushLen(fp))                       return FALSE;
  4099.       hookup(p, q, defaultz);
  4100.       }
  4101.  
  4102.     return TRUE;
  4103.   } /* addElementLen */
  4104.  
  4105.  
  4106. int saveTreeCom (comstrp)
  4107.     char  **comstrp;
  4108.   { /* saveTreeCom */
  4109.     int  ch;
  4110.     boolean  inquote;
  4111.  
  4112.     inquote = FALSE;
  4113.     while ((ch = getc(INFILE)) != EOF && (inquote || ch != ']')) {
  4114.       *(*comstrp)++ = ch;                        /* save character  */
  4115.       if (ch == '[' && ! inquote) {              /* comment; find its end */
  4116.         if ((ch = saveTreeCom(comstrp)) == EOF)  break;
  4117.         *(*comstrp)++ = ch;                      /* add ] */
  4118.         }
  4119.       else if (ch == '\'') inquote = ! inquote;  /* start or end of quote */
  4120.       }
  4121.  
  4122.     return  ch;
  4123.   } /* saveTreeCom */
  4124.  
  4125.  
  4126. boolean processTreeCom (fp, tr)
  4127.     FILE  *fp;
  4128.     tree  *tr;
  4129.   { /* processTreeCom */
  4130.     int   text_started, functor_read, com_open;
  4131.  
  4132.     /*  Accept prefatory "phylip_tree(" or "pseudoNewick("  */
  4133.  
  4134.     functor_read = text_started = 0;
  4135.     (void) fscanf(fp, " p%nhylip_tree(%n", & text_started, & functor_read);
  4136.     if (text_started && ! functor_read) {
  4137.       (void) fscanf(fp, "seudoNewick(%n", & functor_read);
  4138.       if (! functor_read) {
  4139.         printf("Start of tree 'p...' not understood.\n");
  4140.         return FALSE;
  4141.         }
  4142.       }
  4143.  
  4144.     com_open = 0;
  4145.     (void) fscanf(fp, " [%n", & com_open);
  4146.  
  4147.     if (com_open) {                                  /* comment; read it */
  4148.       char  com[1024], *com_end;
  4149.  
  4150.       com_end = com;
  4151.       if (treeFinishCom(fp, & com_end) == EOF) {     /* omits enclosing []s */
  4152.         printf("Missing end of tree comment\n");
  4153.         return FALSE;
  4154.         }
  4155.  
  4156.       (void) readKeyValue(com, likelihood_key, "%lg",
  4157.                                (void *) &(tr->likelihood));
  4158.       (void) readKeyValue(com, opt_level_key,  "%d",
  4159.                                (void *) &(tr->opt_level));
  4160.       (void) readKeyValue(com, smoothed_key,   "%d",
  4161.                                (void *) &(tr->smoothed));
  4162.  
  4163.       if (functor_read) (void) fscanf(fp, " ,");   /* remove trailing comma */
  4164.       }
  4165.  
  4166.     return (functor_read > 0);
  4167.   } /* processTreeCom */
  4168.  
  4169.  
  4170. nodeptr uprootTree (tr, p)
  4171.     tree   *tr;
  4172.     nodeptr p;
  4173.   { /* uprootTree */
  4174.     nodeptr  q, r, s, start;
  4175.     int      n;
  4176.  
  4177.     if (p->tip || p->back) {
  4178.       printf("ERROR: Unable to uproot tree.\n");
  4179.       printf("       Inappropriate node marked for removal.\n");
  4180.       return (nodeptr) NULL;
  4181.       }
  4182.  
  4183.     n = --(tr->nextnode);               /* last internal node added */
  4184.     if (n != tr->mxtips + tr->ntips - 1) {
  4185.       printf("ERROR: Unable to uproot tree.  Inconsistent\n");
  4186.       printf("       number of tips and nodes for rooted tree.\n");
  4187.       return (nodeptr) NULL;
  4188.       }
  4189.  
  4190.     q = p->next->back;                  /* remove p from tree */
  4191.     r = p->next->next->back;
  4192.     hookup(q, r, tr->userlen ? (q->z * r->z) : defaultz);
  4193.  
  4194.     start = (r->tip || (! q->tip)) ? r : r->next->next->back;
  4195.  
  4196.     if (tr->ntips > 2 && p->number != n) {
  4197.       q = tr->nodep[n];            /* transfer last node's conections to p */
  4198.       r = q->next;
  4199.       s = q->next->next;
  4200.       hookup(p,             q->back, q->z);   /* move connections to p */
  4201.       hookup(p->next,       r->back, r->z);
  4202.       hookup(p->next->next, s->back, s->z);
  4203.       if (start->number == q->number) start = start->back->back;
  4204.       q->back = r->back = s->back = (nodeptr) NULL;
  4205.       }
  4206.     else {
  4207.       p->back = p->next->back = p->next->next->back = (nodeptr) NULL;
  4208.       }
  4209.  
  4210.     tr->rooted = FALSE;
  4211.     return  start;
  4212.   } /* uprootTree */
  4213.  
  4214.  
  4215. boolean treeReadLen (fp, tr)
  4216.     FILE  *fp;
  4217.     tree  *tr;
  4218.   { /* treeReadLen */
  4219.     nodeptr  p;
  4220.     int      i, ch;
  4221.     boolean  is_fact;
  4222.  
  4223.     for (i = 1; i <= tr->mxtips; i++) tr->nodep[i]->back = (node *) NULL;
  4224.     tr->start       = tr->nodep[tr->mxtips];
  4225.     tr->ntips       = 0;
  4226.     tr->nextnode    = tr->mxtips + 1;
  4227.     tr->opt_level   = 0;
  4228.     tr->log_f_valid = 0;
  4229.     tr->smoothed    = FALSE;
  4230.     tr->rooted      = FALSE;
  4231.  
  4232.     is_fact = processTreeCom(fp, tr);
  4233.  
  4234.     p = tr->nodep[(tr->nextnode)++];
  4235.     if (! treeNeedCh(fp, '(', "at start of"))       return FALSE;
  4236.     if (! addElementLen(fp, tr, p))                 return FALSE;
  4237.     if (! treeNeedCh(fp, ',', "in"))                return FALSE;
  4238.     if (! addElementLen(fp, tr, p->next))           return FALSE;
  4239.     if (! tr->rooted) {
  4240.       if ((ch = treeGetCh(fp)) == ',') {        /*  An unrooted format */
  4241.         if (! addElementLen(fp, tr, p->next->next)) return FALSE;
  4242.         }
  4243.       else {                                    /*  A rooted format */
  4244.         tr->rooted = TRUE;
  4245.         if (ch != EOF)  (void) ungetc(ch, fp);
  4246.         }
  4247.       }
  4248.     else {
  4249.       p->next->next->back = (nodeptr) NULL;
  4250.       }
  4251.     if (! treeNeedCh(fp, ')', "in"))                return FALSE;
  4252.     (void) treeFlushLabel(fp);
  4253.     if (! treeFlushLen(fp))                         return FALSE;
  4254.     if (is_fact) {
  4255.       if (! treeNeedCh(fp, ')', "at end of"))       return FALSE;
  4256.       if (! treeNeedCh(fp, '.', "at end of"))       return FALSE;
  4257.       }
  4258.     else {
  4259.       if (! treeNeedCh(fp, ';', "at end of"))       return FALSE;
  4260.       }
  4261.  
  4262.     if (tr->rooted) {
  4263.       p->next->next->back = (nodeptr) NULL;
  4264.       tr->start = uprootTree(tr, p->next->next);
  4265.       if (! tr->start)                              return FALSE;
  4266.       }
  4267.     else {
  4268.       tr->start = p->next->next->back;  /* This is start used by treeString */
  4269.       }
  4270.  
  4271.     return  (initrav(tr, tr->start) && initrav(tr, tr->start->back));
  4272.   } /* treeReadLen */
  4273.  
  4274.  
  4275. /*=======================================================================*/
  4276. /*                        Read a tree from a string                      */
  4277. /*=======================================================================*/
  4278.  
  4279.  
  4280. #if Master || Slave
  4281. int str_treeFinishCom (treestrp, strp)
  4282.     char  **treestrp;     /*  tree string pointer */
  4283.     char  **strp;         /*  comment string pointer */
  4284.   { /* str_treeFinishCom */
  4285.     int  ch;
  4286.  
  4287.     while ((ch = *(*treestrp)++) != NULL && ch != ']') {
  4288.       if (strp != NULL) *(*strp)++ = ch;    /* save character  */
  4289.       if (ch == '[') {                      /* nested comment; find its end */
  4290.         if ((ch = str_treeFinishCom(treestrp)) == NULL)  break;
  4291.         if (strp != NULL) *(*strp)++ = ch;  /* save closing ]  */
  4292.         }
  4293.       }
  4294.     if (strp != NULL) **strp = '\0';        /* terminate string  */
  4295.     return  ch;
  4296.   } /* str_treeFinishCom */
  4297.  
  4298.  
  4299. int str_treeGetCh (treestrp)
  4300.     char **treestrp;  /*  tree string pointer */
  4301.     /* get next nonblank, noncomment character */
  4302.   { /* str_treeGetCh */
  4303.     int  ch;
  4304.     
  4305.     while ((ch = *(*treestrp)++) != NULL) {
  4306.       if (whitechar(ch)) ;
  4307.       else if (ch == '[') {                  /* comment; find its end */
  4308.         if ((ch = str_treeFinishCom(treestrp, (char *) NULL)) == NULL)  break;
  4309.         }
  4310.       else  break;
  4311.       }
  4312.  
  4313.     return  ch;
  4314.   } /* str_treeGetCh */
  4315.  
  4316.  
  4317. boolean  str_treeGetLabel (treestrp, lblPtr, maxlen)
  4318.     char **treestrp;  /*  tree string pointer */
  4319.     char  *lblPtr;
  4320.     int    maxlen;
  4321.   { /* str_treeGetLabel */
  4322.     int      ch;
  4323.     boolean  done, quoted, lblfound;
  4324.  
  4325.     if (--maxlen < 0) lblPtr = (char *) NULL;  /* reserves space for '\0' */
  4326.     else if (lblPtr == NULL) maxlen = 0;
  4327.  
  4328.     ch = *(*treestrp)++;
  4329.     done = treeLabelEnd(ch);
  4330.  
  4331.     lblfound = ! done;
  4332.     quoted = (ch == '\'');
  4333.     if (quoted && ! done) {ch = *(*treestrp)++; done = (ch == '\0');}
  4334.  
  4335.     while (! done) {
  4336.       if (quoted) {
  4337.         if (ch == '\'') {ch = *(*treestrp)++; if (ch != '\'') break;}
  4338.         }
  4339.  
  4340.       else if (treeLabelEnd(ch)) break;
  4341.  
  4342.       else if (ch == '_') ch = ' ';  /* unquoted _ goes to space */
  4343.  
  4344.       if (--maxlen >= 0) *lblPtr++ = ch;
  4345.       ch = *(*treestrp)++;
  4346.       if (ch == '\0') break;
  4347.       }
  4348.  
  4349.     (*treestrp)--;
  4350.  
  4351.     if (lblPtr != NULL) *lblPtr = '\0';
  4352.  
  4353.     return lblfound;
  4354.   } /* str_treeGetLabel */
  4355.  
  4356.  
  4357. boolean  str_treeFlushLabel (treestrp)
  4358.     char **treestrp;  /*  tree string pointer */
  4359.   { /* str_treeFlushLabel */
  4360.     return  str_treeGetLabel(treestrp, (char *) NULL, (int) 0);
  4361.   } /* str_treeFlushLabel */
  4362.  
  4363.  
  4364. int  str_treeFindTipName (treestrp, tr)
  4365.     char **treestrp;  /*  tree string pointer */
  4366.     tree  *tr;
  4367.   { /* str_treeFindTipName */
  4368.     nodeptr  q;
  4369.     char    *nameptr, str[nmlngth+2];
  4370.     int      ch, i, n;
  4371.  
  4372.     if (tr->prelabeled) {
  4373.       if (str_treeGetLabel(treestrp, str, nmlngth+2))
  4374.         n = treeFindTipByLabel(str, tr);
  4375.       else
  4376.         n = 0;
  4377.       }
  4378.  
  4379.     else if (tr->ntips < tr->mxtips) {
  4380.       n = tr->ntips + 1;
  4381.       nameptr = tr->nodep[n]->name;
  4382.       if (! str_treeGetLabel(treestrp, nameptr, nmlngth+1)) n = 0;
  4383.       }
  4384.  
  4385.     else {
  4386.       n = 0;
  4387.       }
  4388.  
  4389.     return  n;
  4390.   } /* str_treeFindTipName */
  4391.  
  4392.  
  4393. boolean str_treeProcessLength (treestrp, dptr)
  4394.     char   **treestrp;   /*  tree string ponter */
  4395.     double  *dptr;
  4396.   { /* str_treeProcessLength */
  4397.     int     used;
  4398.  
  4399.     if(! str_treeGetCh(treestrp))  return FALSE;    /*  Skip comments */
  4400.     (*treestrp)--;
  4401.  
  4402.     if (sscanf(*treestrp, "%lf%n", dptr, & used) != 1) {
  4403.       printf("ERROR: str_treeProcessLength: Problem reading branch length\n");
  4404.       printf("%40s\n", *treestrp);
  4405.       *dptr = 0.0;
  4406.       return FALSE;
  4407.       }
  4408.     else {
  4409.       *treestrp += used;
  4410.       }
  4411.  
  4412.     return  TRUE;
  4413.   } /* str_treeProcessLength */
  4414.  
  4415.  
  4416. boolean  str_treeFlushLen (treestrp)
  4417.     char **treestrp;   /*  tree string ponter */
  4418.   { /* str_treeFlushLen */
  4419.     int  ch;
  4420.  
  4421.     if ((ch = str_treeGetCh(treestrp)) == ':')
  4422.       return str_treeProcessLength(treestrp, (double *) NULL);
  4423.     else {
  4424.       (*treestrp)--;
  4425.       return TRUE;
  4426.       }
  4427.   } /* str_treeFlushLen */
  4428.  
  4429.  
  4430. boolean  str_treeNeedCh (treestrp, c1, where)
  4431.     char **treestrp;   /*  tree string pointer */
  4432.     int    c1;
  4433.     char  *where;
  4434.   { /* str_treeNeedCh */
  4435.     int  c2, i;
  4436.  
  4437.     if ((c2 = str_treeGetCh(treestrp)) == c1)  return TRUE;
  4438.  
  4439.     printf("ERROR: Missing '%c' %s tree; ", c1, where);
  4440.     if (c2 == '\0') 
  4441.       printf("end-of-string");
  4442.     else {
  4443.       putchar('"');
  4444.       for (i = 24; i-- && (c2 != '\0'); c2 = *(*treestrp)++)  putchar(c2);
  4445.       putchar('"');
  4446.       }
  4447.  
  4448.     printf(" found instead\n");
  4449.     return FALSE;
  4450.   } /* str_treeNeedCh */
  4451.  
  4452.  
  4453. boolean  str_addElementLen (treestrp, tr, p)
  4454.     char **treestrp;  /*  tree string pointer */
  4455.     tree    *tr;
  4456.     nodeptr  p;
  4457.   { /* str_addElementLen */
  4458.     double   z, branch;
  4459.     nodeptr  q;
  4460.     int      n, ch;
  4461.  
  4462.     if ((ch = str_treeGetCh(treestrp)) == '(') { /*  A new internal node */
  4463.       n = (tr->nextnode)++;
  4464.       if (n > 2*(tr->mxtips) - 2) {
  4465.         if (tr->rooted || n > 2*(tr->mxtips) - 1) {
  4466.           printf("ERROR: too many internal nodes.  Is tree rooted?\n");
  4467.           printf("Deepest splitting should be a trifurcation.\n");
  4468.           return  FALSE;
  4469.           }
  4470.         else {
  4471.           tr->rooted = TRUE;
  4472.           }
  4473.         }
  4474.       q = tr->nodep[n];
  4475.       if (! str_addElementLen(treestrp, tr, q->next))          return FALSE;
  4476.       if (! str_treeNeedCh(treestrp, ',', "in"))               return FALSE;
  4477.       if (! str_addElementLen(treestrp, tr, q->next->next))    return FALSE;
  4478.       if (! str_treeNeedCh(treestrp, ')', "in"))               return FALSE;
  4479.       if (! str_treeFlushLabel(treestrp))                      return FALSE;
  4480.       }
  4481.  
  4482.     else {                           /*  A new tip */
  4483.       n = str_treeFindTipName(treestrp, tr, ch);
  4484.       if (n <= 0) return FALSE;
  4485.       q = tr->nodep[n];
  4486.       if (tr->start->number > n)  tr->start = q;
  4487.       (tr->ntips)++;
  4488.       }                              /* End of tip processing */
  4489.  
  4490.     /*  Master and Slave always use lengths */
  4491.  
  4492.     if (! str_treeNeedCh(treestrp, ':', "in"))                 return FALSE;
  4493.     if (! str_treeProcessLength(treestrp, & branch))           return FALSE;
  4494.     z = exp(-branch / tr->rdta->fracchange);
  4495.     if (z > zmax)  z = zmax;
  4496.     hookup(p, q, z);
  4497.  
  4498.     return  TRUE;
  4499.   } /* str_addElementLen */
  4500.  
  4501.  
  4502. boolean str_processTreeCom(tr, treestrp)
  4503.     tree   *tr;
  4504.     char  **treestrp;
  4505.   { /* str_processTreeCom */
  4506.     char  *com, *com_end;
  4507.     int  text_started, functor_read, com_open;
  4508.  
  4509.     com = *treestrp;
  4510.  
  4511.     functor_read = text_started = 0;
  4512.     sscanf(com, " p%nhylip_tree(%n", & text_started, & functor_read);
  4513.     if (functor_read) {
  4514.       com += functor_read;
  4515.       }
  4516.     else if (text_started) {
  4517.       com += text_started;
  4518.       sscanf(com, "seudoNewick(%n", & functor_read);
  4519.       if (! functor_read) {
  4520.         printf("Start of tree 'p...' not understood.\n");
  4521.         return  FALSE;
  4522.         }
  4523.       else {
  4524.         com += functor_read;
  4525.         }
  4526.       }
  4527.  
  4528.     com_open = 0;
  4529.     sscanf(com, " [%n", & com_open);
  4530.     com += com_open;
  4531.  
  4532.     if (com_open) {                              /* comment; read it */
  4533.     if (!(com_end = strchr(com, ']'))) {
  4534.         printf("Missing end of tree comment.\n");
  4535.         return  FALSE;
  4536.         }
  4537.  
  4538.       *com_end = 0;
  4539.       (void) readKeyValue(com, likelihood_key, "%lg",
  4540.                                (void *) &(tr->likelihood));
  4541.       (void) readKeyValue(com, opt_level_key,  "%d",
  4542.                                (void *) &(tr->opt_level));
  4543.       (void) readKeyValue(com, smoothed_key,   "%d",
  4544.                                (void *) &(tr->smoothed));
  4545.       *com_end = ']';
  4546.       com_end++;
  4547.  
  4548.       if (functor_read) {                          /* remove trailing comma */
  4549.         text_started = 0;
  4550.         sscanf(com_end, " ,%n", & text_started);
  4551.         com_end += text_started;
  4552.         }
  4553.  
  4554.       *treestrp = com_end;
  4555.       }
  4556.  
  4557.     return (functor_read > 0);
  4558.   } /* str_processTreeCom */
  4559.  
  4560.  
  4561. boolean str_treeReadLen (treestr, tr)
  4562.     char  *treestr;
  4563.     tree  *tr;
  4564.     /* read string with representation of tree */
  4565.   { /* str_treeReadLen */
  4566.     nodeptr  p;
  4567.     int  i;
  4568.     boolean  is_fact, found;
  4569.  
  4570.     for (i = 1; i <= tr->mxtips; i++) tr->nodep[i]->back = (node *) NULL;
  4571.     tr->start       = tr->nodep[tr->mxtips];
  4572.     tr->ntips       = 0;
  4573.     tr->nextnode    = tr->mxtips + 1;
  4574.     tr->opt_level   = 0;
  4575.     tr->log_f_valid = 0;
  4576.     tr->smoothed    = Master;
  4577.     tr->rooted      = FALSE;
  4578.  
  4579.     is_fact = str_processTreeCom(tr, & treestr);
  4580.  
  4581.     p = tr->nodep[(tr->nextnode)++];
  4582.     if (! str_treeNeedCh(& treestr, '(', "at start of"))       return FALSE;
  4583.     if (! str_addElementLen(& treestr, tr, p))                 return FALSE;
  4584.     if (! str_treeNeedCh(& treestr, ',', "in"))                return FALSE;
  4585.     if (! str_addElementLen(& treestr, tr, p->next))           return FALSE;
  4586.     if (! tr->rooted) {
  4587.       if (str_treeGetCh(& treestr) == ',') {        /*  An unrooted format */
  4588.         if (! str_addElementLen(& treestr, tr, p->next->next)) return FALSE;
  4589.         }
  4590.       else {                                       /*  A rooted format */
  4591.         p->next->next->back = (nodeptr) NULL;
  4592.         tr->rooted = TRUE;
  4593.         treestr--;
  4594.         }
  4595.       }
  4596.     if (! str_treeNeedCh(& treestr, ')', "in"))                return FALSE;
  4597.     if (! str_treeFlushLabel(& treestr))                       return FALSE;
  4598.     if (! str_treeFlushLen(& treestr))                         return FALSE;
  4599.     if (is_fact) {
  4600.       if (! str_treeNeedCh(& treestr, ')', "at end of"))       return FALSE;
  4601.       if (! str_treeNeedCh(& treestr, '.', "at end of"))       return FALSE;
  4602.       }
  4603.     else {
  4604.       if (! str_treeNeedCh(& treestr, ';', "at end of"))       return FALSE;
  4605.       }
  4606.  
  4607.     if (tr->rooted)  if (! uprootTree(tr, p->next->next))     return FALSE;
  4608.     tr->start = p->next->next->back;  /* This is start used by treeString */
  4609.  
  4610.     return  (initrav(tr, tr->start) && initrav(tr, tr->start->back));
  4611.   } /* str_treeReadLen */
  4612. #endif
  4613.  
  4614.  
  4615. boolean treeEvaluate (tr, bt)         /* Evaluate a user tree */
  4616.     tree     *tr;
  4617.     bestlist *bt;
  4618.   { /* treeEvaluate */
  4619.  
  4620.     if (Slave || ! tr->userlen) {
  4621.       if (! smoothTree(tr, 4 * smoothings)) return FALSE;
  4622.       }
  4623.  
  4624.     if (evaluate(tr, tr->start) == badEval)  return FALSE;
  4625.  
  4626. #   if ! Slave
  4627.       (void) saveBestTree(bt, tr);
  4628. #   endif
  4629.     return TRUE;
  4630.   } /* treeEvaluate */
  4631.  
  4632.  
  4633. #if Master || Slave
  4634. FILE *freopen_pid (filenm, mode, stream)
  4635.     char *filenm, *mode;
  4636.     FILE *stream;
  4637.   { /* freopen_pid */
  4638.     char scr[512];
  4639.  
  4640.     (void) sprintf(scr, "%s.%d", filenm, getpid());
  4641.     return  freopen(scr, mode, stream);
  4642.   } /* freopen_pid */
  4643. #endif
  4644.  
  4645.  
  4646. boolean  showBestTrees (bt, tr, adef, treefile)
  4647.     bestlist *bt;
  4648.     tree     *tr;
  4649.     analdef  *adef;
  4650.     FILE     *treefile;
  4651.   { /* showBestTrees */
  4652.     int     rank;
  4653.  
  4654.     for (rank = 1; rank <= bt->nvalid; rank++) {
  4655.       if (rank > 1) {
  4656.         if (rank != recallBestTree(bt, rank, tr))  break;
  4657.         }
  4658.       if (evaluate(tr, tr->start) == badEval) return FALSE;
  4659.       if (tr->outgrnode->back)  tr->start = tr->outgrnode;
  4660.       printTree(tr, adef);
  4661.       summarize(tr);
  4662.       if (treefile)  treeOut(treefile, tr, adef->trout);
  4663.       }
  4664.  
  4665.     return TRUE;
  4666.   } /* showBestTrees */
  4667.  
  4668.  
  4669. boolean cmpBestTrees (bt, tr)
  4670.     bestlist *bt;
  4671.     tree     *tr;
  4672.   { /* cmpBestTrees */
  4673.     double  sum, sum2, sd, temp, wtemp, bestscore;
  4674.     double *log_f0, *log_f0_ptr;      /* Save a copy of best log_f */
  4675.     double *log_f_ptr;
  4676.     int     i, j, num, besttips;
  4677.  
  4678.     num = bt->nvalid;
  4679.     if ((num <= 1) || (tr->cdta->wgtsum <= 1)) return TRUE;
  4680.  
  4681.     if (! (log_f0 = (double *) Malloc(sizeof(double) * tr->cdta->endsite))) {
  4682.       printf("ERROR: cmpBestTrees unable to obtain space for log_f0\n");
  4683.       return FALSE;
  4684.       }
  4685.  
  4686.     printf("Tree      Ln L        Diff Ln L       Its S.D.");
  4687.     printf("   Significantly worse?\n\n");
  4688.  
  4689.     for (i = 1; i <= num; i++) {
  4690.       if (i != recallBestTree(bt, i, tr))  break;
  4691.       if (! (tr->log_f_valid))  {
  4692.         if (evaluate(tr, tr->start) == badEval) return FALSE;
  4693.         }
  4694.  
  4695.       printf("%3d%14.5f", i, tr->likelihood);
  4696.       if (i == 1) {
  4697.         printf("  <------ best\n");
  4698.         besttips = tr->ntips;
  4699.         bestscore = tr->likelihood;
  4700.         log_f0_ptr = log_f0;
  4701.         log_f_ptr  = tr->log_f;
  4702.         for (j = 0; j < tr->cdta->endsite; j++)  *log_f0_ptr++ = *log_f_ptr++;
  4703.         }
  4704.       else if (tr->ntips != besttips)
  4705.         printf("  (different number of species)\n");
  4706.       else {
  4707.         sum = sum2 = 0.0;
  4708.         log_f0_ptr = log_f0;
  4709.         log_f_ptr  = tr->log_f;
  4710.         for (j = 0; j < tr->cdta->endsite; j++) {
  4711.           temp  = *log_f0_ptr++ - *log_f_ptr++;
  4712.           wtemp = tr->cdta->aliaswgt[j] * temp;
  4713.           sum  += wtemp;
  4714.           sum2 += wtemp * temp;
  4715.           }
  4716.         sd = sqrt( tr->cdta->wgtsum * (sum2 - sum*sum / tr->cdta->wgtsum)
  4717.                                     / (tr->cdta->wgtsum - 1) );
  4718.         printf("%14.5f%14.4f", tr->likelihood - bestscore, sd);
  4719.         printf("           %s\n", (sum > 1.95996 * sd) ? "Yes" : " No");
  4720.         }
  4721.       }
  4722.  
  4723.     Free(log_f0);
  4724.     printf("\n\n");
  4725.  
  4726.     return TRUE;
  4727.   } /* cmpBestTrees */
  4728.  
  4729.  
  4730. boolean  makeUserTree (tr, bt, adef)
  4731.     tree     *tr;
  4732.     bestlist *bt;
  4733.     analdef  *adef;
  4734.   { /* makeUserTree */
  4735.     char   filename[128];
  4736.     FILE  *treefile;
  4737.     int    nusertrees, which;
  4738.  
  4739.     nusertrees = adef->numutrees;
  4740.  
  4741.     printf("User-defined %s:\n\n", (nusertrees == 1) ? "tree" : "trees");
  4742.  
  4743.     treefile = adef->trout ? fopen_pid("treefile", "w", filename) : (FILE *) NULL;
  4744.  
  4745.     for (which = 1; which <= nusertrees; which++) {
  4746.       if (! treeReadLen(INFILE, tr)) return FALSE;
  4747.       if (! treeEvaluate(tr, bt))    return FALSE;
  4748.       if (tr->global <= 0) {
  4749.         if (tr->outgrnode->back)  tr->start = tr->outgrnode;
  4750.         printTree(tr, adef);
  4751.         summarize(tr);
  4752.         if (treefile)  treeOut(treefile, tr, adef->trout);
  4753.         }
  4754.       else {
  4755.         printf("%6d:  Ln Likelihood =%14.5f\n", which, tr->likelihood);
  4756.         }
  4757.       }
  4758.  
  4759.     if (tr->global > 0) {
  4760.       putchar('\n');
  4761.       if (! recallBestTree(bt, 1, tr))  return FALSE;
  4762.       printf("      Ln Likelihood =%14.5f\n", tr->likelihood);
  4763.       if (! optimize(tr, tr->global, bt))  return FALSE;
  4764.       if (tr->outgrnode->back)  tr->start = tr->outgrnode;
  4765.       printTree(tr, adef);
  4766.       summarize(tr);
  4767.       if (treefile)  treeOut(treefile, tr, adef->trout);
  4768.       }
  4769.  
  4770.     if (treefile) {
  4771.       (void) fclose(treefile);
  4772.       printf("Tree also written to %s\n", filename);
  4773.       }
  4774.  
  4775.     putchar('\n');
  4776.  
  4777.     (void) cmpBestTrees(bt, tr);
  4778.     return TRUE;
  4779.   } /* makeUserTree */
  4780.  
  4781.  
  4782. #if Slave
  4783. boolean slaveTreeEvaluate (tr, bt)
  4784.     tree     *tr;
  4785.     bestlist *bt;
  4786.   { /* slaveTreeEvaluate */
  4787.     boolean done;
  4788.  
  4789.     do {
  4790.        requestForWork();
  4791.        if (! receiveTree(& comm_master, tr))        return FALSE;
  4792.        done = tr->likelihood > 0.0;
  4793.        if (! done) {
  4794.          if (! treeEvaluate(tr, bt))                return FALSE;
  4795.          if (! sendTree(& comm_master, tr))         return FALSE;
  4796.          }
  4797.        } while (! done);
  4798.  
  4799.     return TRUE;
  4800.   } /* slaveTreeEvaluate */
  4801. #endif
  4802.  
  4803.  
  4804. double randum (seed)
  4805.     long  *seed;
  4806.     /* random number generator, modified to use 12 bit chunks */
  4807.   { /* randum */
  4808.     long  sum, mult0, mult1, seed0, seed1, seed2, newseed0, newseed1, newseed2;
  4809.  
  4810.     mult0 = 1549;
  4811.     seed0 = *seed & 4095;
  4812.     sum  = mult0 * seed0;
  4813.     newseed0 = sum & 4095;
  4814.     sum >>= 12;
  4815.     seed1 = (*seed >> 12) & 4095;
  4816.     mult1 =  406;
  4817.     sum += mult0 * seed1 + mult1 * seed0;
  4818.     newseed1 = sum & 4095;
  4819.     sum >>= 12;
  4820.     seed2 = (*seed >> 24) & 255;
  4821.     sum += mult0 * seed2 + mult1 * seed1;
  4822.     newseed2 = sum & 255;
  4823.  
  4824.     *seed = newseed2 << 24 | newseed1 << 12 | newseed0;
  4825.     return  0.00390625 * (newseed2
  4826.                           + 0.000244140625 * (newseed1
  4827.                                               + 0.000244140625 * newseed0));
  4828.   } /* randum */
  4829.  
  4830.  
  4831. boolean makeDenovoTree (tr, bt, adef)
  4832.     tree     *tr;
  4833.     bestlist *bt;
  4834.     analdef  *adef;
  4835.   { /* makeDenovoTree */
  4836.     char   filename[128];
  4837.     FILE  *treefile;
  4838.     nodeptr  p;
  4839.     int  *enterorder;      /*  random entry order */
  4840.     int  i, j, k, nextsp, newsp, maxtrav, tested;
  4841.  
  4842.     double randum();
  4843.  
  4844.  
  4845.     enterorder = (int *) Malloc(sizeof(int) * (tr->mxtips + 1));
  4846.     if (! enterorder) {
  4847.        fprintf(stderr, "makeDenovoTree: Malloc failure for enterorder\n");
  4848.        return 0;
  4849.        }
  4850.  
  4851.     if (adef->restart) {
  4852.       printf("Restarting from tree with the following sequences:\n");
  4853.       tr->userlen = TRUE;
  4854.       if (! treeReadLen(INFILE, tr))          return FALSE;
  4855.       if (! smoothTree(tr, smoothings))       return FALSE;
  4856.       if (evaluate(tr, tr->start) == badEval) return FALSE;
  4857.       if (saveBestTree(bt, tr) < 1)           return FALSE;
  4858.  
  4859.       for (i = 1, j = tr->ntips; i <= tr->mxtips; i++) { /* find loose tips */
  4860.         if (! tr->nodep[i]->back) {
  4861.           enterorder[++j] = i;
  4862.           }
  4863.         else {
  4864.           printf("   %s\n", tr->nodep[i]->name);
  4865.  
  4866. #         if Master
  4867.             if (i>3) REPORT_ADD_SPECS;
  4868. #         endif
  4869.           }
  4870.         }
  4871.       putchar('\n');
  4872.       }
  4873.  
  4874.     else {                                           /* start from scratch */
  4875.       tr->ntips = 0;
  4876.       for (i = 1; i <= tr->mxtips; i++) enterorder[i] = i;
  4877.       }
  4878.  
  4879.     if (adef->jumble) for (i = tr->ntips + 1; i <= tr->mxtips; i++) {
  4880.       j = randum(&(adef->jumble))*(tr->mxtips - tr->ntips) + tr->ntips + 1;
  4881.       k = enterorder[j];
  4882.       enterorder[j] = enterorder[i];
  4883.       enterorder[i] = k;
  4884.       }
  4885.  
  4886.     bt->numtrees = 1;
  4887.     if (tr->ntips < tr->mxtips)  printf("Adding species:\n");
  4888.  
  4889.     if (tr->ntips == 0) {
  4890.       for (i = 1; i <= 3; i++) {
  4891.         printf("   %s\n", tr->nodep[enterorder[i]]->name);
  4892.         }
  4893.       tr->nextnode = tr->mxtips + 1;
  4894.       if (! buildSimpleTree(tr, enterorder[1], enterorder[2], enterorder[3]))
  4895.         return FALSE;
  4896.       }
  4897.  
  4898.     while (tr->ntips < tr->mxtips || tr->opt_level < tr->global) {
  4899.       maxtrav = (tr->ntips == tr->mxtips) ? tr->global : tr->partswap;
  4900.       if (maxtrav > tr->ntips - 3)  maxtrav = tr->ntips - 3;
  4901.  
  4902.       if (tr->opt_level >= maxtrav) {
  4903.         nextsp = ++(tr->ntips);
  4904.         newsp = enterorder[nextsp];
  4905.         p = tr->nodep[newsp];
  4906.         printf("   %s\n", p->name);
  4907.  
  4908. #       if Master
  4909.           if (nextsp % DNAML_STEP_TIME_COUNT == 1) {
  4910.             REPORT_STEP_TIME;
  4911.             }
  4912.           REPORT_ADD_SPECS;
  4913. #       endif
  4914.  
  4915.         (void) buildNewTip(tr, p);
  4916.  
  4917.         resetBestTree(bt);
  4918.         cacheZ(tr);
  4919.         tested = addTraverse(tr, p->back, findAnyTip(tr->start)->back,
  4920.                              1, tr->ntips - 2, bt, adef->qadd);
  4921.         if (tested == badRear) return FALSE;
  4922.         bt->numtrees += tested;
  4923.  
  4924. #       if Master
  4925.           getReturnedTrees(tr, bt, tested);
  4926. #       endif
  4927.  
  4928.         printf("      Tested %d alternative trees\n", tested);
  4929.  
  4930.         (void) recallBestTree(bt, 1, tr);
  4931.         if (! tr->smoothed) {
  4932.           if (! smoothTree(tr, smoothings))        return FALSE;
  4933.           if (evaluate(tr, tr->start) == badEval)  return FALSE;
  4934.           (void) saveBestTree(bt, tr);
  4935.           }
  4936.  
  4937.         if (tr->ntips == 4)  tr->opt_level = 1;  /* All 4 taxon trees done */
  4938.         maxtrav = (tr->ntips == tr->mxtips) ? tr->global : tr->partswap;
  4939.         if (maxtrav > tr->ntips - 3)  maxtrav = tr->ntips - 3;
  4940.         }
  4941.  
  4942.       printf("      Ln Likelihood =%14.5f\n", tr->likelihood);
  4943.       if (! optimize(tr, maxtrav, bt)) return FALSE;
  4944.       }
  4945.  
  4946.     printf("\nExamined %d %s\n", bt->numtrees,
  4947.                                  bt->numtrees != 1 ? "trees" : "tree");
  4948.  
  4949.     treefile = adef->trout ? fopen_pid("treefile", "w", filename) : (FILE *) NULL;
  4950.     (void) showBestTrees(bt, tr, adef, treefile);
  4951.     if (treefile) {
  4952.       (void) fclose(treefile);
  4953.       printf("Tree also written to %s\n\n", filename);
  4954.       }
  4955.  
  4956.     (void) cmpBestTrees(bt, tr);
  4957.  
  4958. #   if DeleteCheckpointFile
  4959.       unlink_pid(checkpointname);
  4960. #   endif
  4961.  
  4962.     Free(enterorder);
  4963.  
  4964.     return TRUE;
  4965.   } /* makeDenovoTree */
  4966.  
  4967. /*==========================================================================*/
  4968. /*                             "main" routine                               */
  4969. /*==========================================================================*/
  4970.  
  4971. #if defined(MACINTOSH) || defined(MSDOS)
  4972. int RealMain(int argc, char** argv)
  4973. #else
  4974. #if Sequential
  4975.   main ()
  4976. #else
  4977.   slave ()
  4978. #endif
  4979. #endif
  4980.   { /* DNA Maximum Likelihood */
  4981. #   if Master
  4982.       int starttime, inputtime, endtime;
  4983. #   endif
  4984.  
  4985. #   if Sequential
  4986.         /* dgg addition */
  4987.       long starttime, inputtime, endtime;
  4988.       double millisecs;
  4989. #   endif
  4990.  
  4991. #   if Master || Slave
  4992.       int my_id, nprocs, type, from, sz;
  4993.       char *msg;
  4994. #   endif
  4995.  
  4996.     analdef      *adef;
  4997.     rawDNA       *rdta;
  4998.     crunchedDNA  *cdta;
  4999.     tree         *tr;    /*  current tree */
  5000.     bestlist     *bt;    /*  topology of best found tree */
  5001.  
  5002.  
  5003. #   if Debug
  5004.       {
  5005.         char debugfilename[128];
  5006.         debug = fopen_pid("dnaml_debug", "w", debugfilename);
  5007.         }
  5008. #   endif
  5009.  
  5010. #   if Master
  5011.       starttime = p4_clock();
  5012.       nprocs = p4_num_total_slaves();
  5013.  
  5014.       if ((OUTFILE = freopen_pid("master.out", "w", stdout)) == NULL) {
  5015.         fprintf(stderr, "Could not open output file\n");
  5016.         exit(1);
  5017.         }
  5018.  
  5019.       /* Receive input file name from host */
  5020.       type = DNAML_FILE_NAME;
  5021.       from = DNAML_HOST_ID;
  5022.       msg  = NULL;
  5023.       p4_recv(& type, & from, & msg, & sz);
  5024.       if ((INFILE = fopen(msg, "r")) == NULL) {
  5025.         fprintf(stderr, "master could not open input file %s\n", msg);
  5026.         exit(1);
  5027.         }
  5028.       p4_msg_free(msg);
  5029.  
  5030.       open_link(& comm_slave);
  5031. #   endif
  5032.  
  5033. #   if Sequential
  5034.       starttime = clock();
  5035. #ifdef NoSTDIN
  5036.             {
  5037.             char *fname = "infile";
  5038.       if ((INFILE = fopen(fname, "r")) == NULL) 
  5039.       do {
  5040.           char filename[512];
  5041.           filename[0]= 0;
  5042.         fprintf(stderr, "Could not open input file '%s'\n", fname);
  5043.         fprintf(stderr, "Enter path-name for data file: ");
  5044.         scanf( "%s", filename);
  5045.         fname= filename;
  5046.            INFILE = fopen(fname, "r");
  5047.       } while (!INFILE && *fname > ' ');
  5048.       if (INFILE == NULL) {
  5049.              fprintf(stderr, "Could not open input file '%s'\n", fname);
  5050.             exit(1);
  5051.           }
  5052.       }
  5053. #endif
  5054. #        endif
  5055.  
  5056. #  if DNAML_STEP
  5057.       begin_step_time = starttime;
  5058. #  endif
  5059.  
  5060. #   if Slave
  5061.       my_id = p4_get_my_id();
  5062.       nprocs = p4_num_total_slaves();
  5063.  
  5064.       /* Receive input file name from host */
  5065.       type = DNAML_FILE_NAME;
  5066.       from = DNAML_HOST_ID;
  5067.       msg  = NULL;
  5068.       p4_recv(& type, & from, & msg, & sz);
  5069.       if ((INFILE = fopen(msg, "r")) == NULL) {
  5070.         fprintf(stderr, "slave could not open input file %s\n",msg);
  5071.         exit(1);
  5072.         }
  5073.       p4_msg_free(msg);
  5074.  
  5075. #     ifdef P4DEBUG
  5076.         if ((OUTFILE = freopen_pid("slave.out", "w", stdout)) == NULL) {
  5077.           fprintf(stderr, "Could not open output file\n");
  5078.           exit(1);
  5079.           }
  5080. #     else
  5081.         if ((OUTFILE = freopen("/dev/null", "w", stdout)) == NULL) {
  5082.           fprintf(stderr, "Could not open output file\n");
  5083.           exit(1);
  5084.           }
  5085. #     endif
  5086.  
  5087.       open_link(& comm_master);
  5088. #   endif
  5089.  
  5090.  
  5091. /*  Get data structure memory  */
  5092.  
  5093.     if (! (adef = (analdef *) Malloc(sizeof(analdef)))) {
  5094.       printf("ERROR: Unable to get memory for analysis definition\n\n");
  5095.       return 1;
  5096.       }
  5097.  
  5098.     if (! (rdta = (rawDNA *) Malloc(sizeof(rawDNA)))) {
  5099.       printf("ERROR: Unable to get memory for raw DNA\n\n");
  5100.       return 1;
  5101.       }
  5102.  
  5103.     if (! (cdta = (crunchedDNA *) Malloc(sizeof(crunchedDNA)))) {
  5104.       printf("ERROR: Unable to get memory for crunched DNA\n\n");
  5105.       return 1;
  5106.       }
  5107.  
  5108.     if ((tr = (tree *)     Malloc(sizeof(tree))) &&
  5109.         (bt = (bestlist *) Malloc(sizeof(bestlist)))) ;
  5110.     else {
  5111.       printf("ERROR: Unable to get memory for trees\n\n");
  5112.       return 1;
  5113.       }
  5114.     bt->ninit = 0;
  5115.  
  5116.     if (! getinput(adef, rdta, cdta, tr))                            return 1;
  5117.  
  5118. #   if Master
  5119.       inputtime = p4_clock();
  5120.       printf("Input time %d milliseconds\n", inputtime - starttime);
  5121.       REPORT_STEP_TIME;
  5122. #   endif
  5123.  
  5124. #   if Slave
  5125.       (void) fclose(INFILE);
  5126. #   endif
  5127.  
  5128. #   if Sequential
  5129.       inputtime = clock();
  5130.       millisecs= TIMETOSECS(inputtime - starttime);
  5131.       printf("Input time %6.4f seconds\n",millisecs);
  5132. #   ifdef NoSTDIN
  5133.       (void) fclose(INFILE);
  5134. #   endif
  5135. #   endif
  5136.  
  5137.  
  5138. /*  The material below would be a loop over jumbles and/or boots */
  5139.  
  5140.     if (! makeweights(adef, rdta, cdta))                             return 1;
  5141.     if (! makevalues(rdta, cdta))                                    return 1;
  5142.     if (adef->empf && ! empiricalfreqs(rdta, cdta))                  return 1;
  5143.     reportfreqs(adef, rdta);
  5144.     if (! linkdata2tree(rdta, cdta, tr))                             return 1;
  5145.  
  5146.     if (! linkxarray(3, 3, cdta->endsite, & freextip, & usedxtip))   return 1;
  5147.     if (! setupnodex(tr))                                            return 1;
  5148.  
  5149. #   if Slave
  5150.       if (! slaveTreeEvaluate(tr, bt))                               return 1;
  5151. #   else
  5152.       if (! initBestTree(bt, adef->nkeep, tr->mxtips, tr->cdta->endsite))    return 1;
  5153.       if (! adef->usertree) {
  5154.         if (! makeDenovoTree(tr, bt, adef))                          return 1;
  5155.         }
  5156.       else {
  5157.         if (! makeUserTree(tr, bt, adef))                            return 1;
  5158.         }
  5159.       if (! freeBestTree(bt))                                        return 1;
  5160. #   endif
  5161.  
  5162. /*  Endpoint for jumble and/or boot loop */
  5163.  
  5164. #   if Master
  5165.       tr->likelihood = 1.0;             /* terminate slaves */
  5166.       (void) sendTree(& comm_slave, tr);
  5167. #   endif
  5168.  
  5169.     freeTree(tr);
  5170.  
  5171. #   if Master
  5172.       close_link(& comm_slave);
  5173.       (void) fclose(INFILE);
  5174.  
  5175.       REPORT_STEP_TIME;
  5176.       endtime = p4_clock();
  5177.       printf("Execution time %d milliseconds\n", endtime - inputtime);
  5178.       (void) fclose(OUTFILE);
  5179. #   endif
  5180.  
  5181. #   if Slave
  5182.       close_link(& comm_master);
  5183.       (void) fclose(OUTFILE);
  5184. #   endif
  5185.  
  5186. #   if Sequential
  5187.         /* dgg addition */
  5188.       endtime = clock();
  5189.       millisecs= TIMETOSECS(endtime - inputtime);
  5190.       printf("Execution time %6.4f seconds\n", millisecs);
  5191. #   endif
  5192.  
  5193. #   if Debug
  5194.       (void) fclose(debug);
  5195. #   endif
  5196.  
  5197. #   if Master || Slave
  5198.       p4_send(DNAML_DONE, DNAML_HOST_ID, NULL, 0);
  5199. #   else
  5200.       return 0;
  5201. #   endif
  5202.   } /* DNA Maximum Likelihood */
  5203.